{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Expected models on binned contact maps\n",
    "\n",
    "\n",
    "* For intrachromosomal arm regions: P(s) by diagonal\n",
    "* For interchromosomal regions: Average contact frequency by block\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:10:52.840179Z",
     "start_time": "2018-05-08T02:10:51.787737Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import seaborn as sns\n",
    "import multiprocess as mp\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import bioframe\n",
    "import cooltools\n",
    "import cooler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:10:53.193847Z",
     "start_time": "2018-05-08T02:10:52.845583Z"
    }
   },
   "outputs": [],
   "source": [
    "mm9 = bioframe.fetch_chromsizes('mm9')\n",
    "chromsizes = bioframe.fetch_chromsizes('mm9')\n",
    "chromosomes = list(chromsizes.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:10:53.220550Z",
     "start_time": "2018-05-08T02:10:53.196035Z"
    }
   },
   "outputs": [],
   "source": [
    "conditions = ['WT', 'T', 'dN']\n",
    "binsize = 100000\n",
    "\n",
    "cooler_paths = {    \n",
    "    'WT' : f'data/UNTR.{binsize}.cool',\n",
    "    'T'  : f'data/TAM.{binsize}.cool',\n",
    "    'dN' : f'data/NIPBL.{binsize}.cool',\n",
    "}\n",
    "long_names = {\n",
    "    'WT': 'Wildtype',\n",
    "    'T' : 'TAM',\n",
    "    'dN': 'NipblKO',\n",
    "}\n",
    "pal = sns.color_palette('colorblind')\n",
    "colors = {\n",
    "    'WT': pal[0],\n",
    "    'T' : '#333333',\n",
    "    'dN': pal[2],\n",
    "}\n",
    "\n",
    "clrs = {\n",
    "    cond: cooler.Cooler(cooler_paths[cond]) for cond in conditions\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:33.313532Z",
     "start_time": "2018-05-08T02:10:53.222362Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WT cis\n",
      "WT trans\n",
      "T cis\n",
      "T trans\n",
      "dN cis\n",
      "dN trans\n"
     ]
    }
   ],
   "source": [
    "# this cell takes a long time to run\n",
    "from cooltools.expected import diagsum, blocksum_pairwise\n",
    "supports = [(chrom, 0, chromsizes[chrom]) for chrom in chromosomes]\n",
    "\n",
    "\n",
    "cis_exp = {}\n",
    "trs_exp = {}\n",
    "\n",
    "with mp.Pool() as pool:\n",
    "    for cond in conditions:\n",
    "        \n",
    "        print(cond, 'cis')\n",
    "        tables = diagsum(\n",
    "            clrs[cond], \n",
    "            supports, \n",
    "            transforms={\n",
    "                'balanced': lambda p: p['count'] * p['weight1'] * p['weight2'],\n",
    "            },\n",
    "            chunksize=10000000,\n",
    "            ignore_diags=2, \n",
    "            map=pool.map)\n",
    "        \n",
    "        cis_exp[cond] = pd.concat(\n",
    "            [tables[support] for support in supports], \n",
    "            keys=[support[0] for support in supports], \n",
    "            names=['chrom'])\n",
    "        cis_exp[cond]['balanced.avg'] = cis_exp[cond]['balanced.sum'] / cis_exp[cond]['n_valid']\n",
    "        \n",
    "        cis_exp[cond].to_csv(f'data/{long_names[cond]}.{binsize//1000}kb.expected.cis.tsv', sep='\\t')\n",
    "        \n",
    "        print(cond, 'trans')\n",
    "        records = blocksum_pairwise(\n",
    "            clrs[cond],                \n",
    "            supports, \n",
    "            transforms={\n",
    "                'balanced': lambda p: p['count'] * p['weight1'] * p['weight2'],\n",
    "            },\n",
    "            chunksize=10000000,\n",
    "            map=pool.map)\n",
    "        \n",
    "        trs_exp[cond] = pd.DataFrame(\n",
    "            [{'chrom1': s1[0], 'chrom2': s2[0], **rec} for (s1, s2), rec in records.items()], \n",
    "            columns=['chrom1', 'chrom2', 'n_valid', 'count.sum', 'balanced.sum'])\n",
    "        \n",
    "        trs_exp[cond].to_csv(f'data/{long_names[cond]}.{binsize//1000}kb.expected.trans.tsv', sep='\\t', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:55.074487Z",
     "start_time": "2018-05-08T02:12:55.061941Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>n_valid</th>\n",
       "      <th>count.sum</th>\n",
       "      <th>balanced.sum</th>\n",
       "      <th>balanced.avg</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>chrom</th>\n",
       "      <th>diag</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">chr1</th>\n",
       "      <th>0</th>\n",
       "      <td>18900</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>18652</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>18638</td>\n",
       "      <td>132754.0</td>\n",
       "      <td>301.561046</td>\n",
       "      <td>0.016180</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>18634</td>\n",
       "      <td>92489.0</td>\n",
       "      <td>208.895766</td>\n",
       "      <td>0.011210</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>18623</td>\n",
       "      <td>70473.0</td>\n",
       "      <td>159.212659</td>\n",
       "      <td>0.008549</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            n_valid  count.sum  balanced.sum  balanced.avg\n",
       "chrom diag                                                \n",
       "chr1  0       18900        NaN           NaN           NaN\n",
       "      1       18652        NaN           NaN           NaN\n",
       "      2       18638   132754.0    301.561046      0.016180\n",
       "      3       18634    92489.0    208.895766      0.011210\n",
       "      4       18623    70473.0    159.212659      0.008549"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cis_exp['WT'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:55.415342Z",
     "start_time": "2018-05-08T02:12:55.404914Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>chrom1</th>\n",
       "      <th>chrom2</th>\n",
       "      <th>n_valid</th>\n",
       "      <th>count.sum</th>\n",
       "      <th>balanced.sum</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>chr1</td>\n",
       "      <td>chr2</td>\n",
       "      <td>357714000</td>\n",
       "      <td>147693.0</td>\n",
       "      <td>338.795457</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>chr1</td>\n",
       "      <td>chr3</td>\n",
       "      <td>314220340</td>\n",
       "      <td>127220.0</td>\n",
       "      <td>309.070464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>chr1</td>\n",
       "      <td>chr4</td>\n",
       "      <td>306122580</td>\n",
       "      <td>123733.0</td>\n",
       "      <td>290.960206</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>chr1</td>\n",
       "      <td>chr5</td>\n",
       "      <td>299987240</td>\n",
       "      <td>112388.0</td>\n",
       "      <td>264.008169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>chr1</td>\n",
       "      <td>chr6</td>\n",
       "      <td>294360620</td>\n",
       "      <td>119053.0</td>\n",
       "      <td>280.474475</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  chrom1 chrom2    n_valid  count.sum  balanced.sum\n",
       "0   chr1   chr2  357714000   147693.0    338.795457\n",
       "1   chr1   chr3  314220340   127220.0    309.070464\n",
       "2   chr1   chr4  306122580   123733.0    290.960206\n",
       "3   chr1   chr5  299987240   112388.0    264.008169\n",
       "4   chr1   chr6  294360620   119053.0    280.474475"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trs_exp['WT'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:56.634905Z",
     "start_time": "2018-05-08T02:12:56.612745Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cis</th>\n",
       "      <th>trans</th>\n",
       "      <th>total</th>\n",
       "      <th>cis:trans</th>\n",
       "      <th>cis:total</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>NipblKO</th>\n",
       "      <td>60644854</td>\n",
       "      <td>29374071</td>\n",
       "      <td>90018925</td>\n",
       "      <td>2.064571</td>\n",
       "      <td>0.673690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>TAM</th>\n",
       "      <td>46889387</td>\n",
       "      <td>19725414</td>\n",
       "      <td>66614801</td>\n",
       "      <td>2.377105</td>\n",
       "      <td>0.703888</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Wildtype</th>\n",
       "      <td>37717628</td>\n",
       "      <td>14296684</td>\n",
       "      <td>52014312</td>\n",
       "      <td>2.638208</td>\n",
       "      <td>0.725139</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               cis     trans     total  cis:trans  cis:total\n",
       "NipblKO   60644854  29374071  90018925   2.064571   0.673690\n",
       "TAM       46889387  19725414  66614801   2.377105   0.703888\n",
       "Wildtype  37717628  14296684  52014312   2.638208   0.725139"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stats = {}\n",
    "for cond in conditions:    \n",
    "    n_cis = int(cis_exp[cond]['count.sum'].sum())\n",
    "    n_trs = int(trs_exp[cond]['count.sum'].sum())\n",
    "    stats[long_names[cond]] = {\n",
    "        'cis': n_cis,\n",
    "        'trans': n_trs,\n",
    "        'total': n_cis + n_trs,\n",
    "        'cis:trans': n_cis / n_trs,\n",
    "        'cis:total': n_cis / (n_cis + n_trs)\n",
    "    }\n",
    "pd.DataFrame.from_dict(stats, orient='index')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:57.412670Z",
     "start_time": "2018-05-08T02:12:57.331551Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/nezar/miniconda3/envs/py36/lib/python3.6/site-packages/ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in true_divide\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "sums = {}\n",
    "n_valid = {}\n",
    "scalings = {}\n",
    "for cond in conditions:\n",
    "    grouped = cis_exp[cond].groupby('diag')\n",
    "    n_valid[cond] = grouped['n_valid'].sum().values\n",
    "    sums[cond] = grouped['balanced.sum'].sum().values\n",
    "    scalings[cond] = (sums[cond] / n_valid[cond])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:12:58.092032Z",
     "start_time": "2018-05-08T02:12:58.082970Z"
    }
   },
   "outputs": [],
   "source": [
    "from cooltools.lib import numutils\n",
    "\n",
    "def coarsen_geometric(sums, counts, n_bins=100):\n",
    "    \"\"\"Re-bin the expected sums into logarithmically growing bins.\n",
    "    \n",
    "    \"\"\"\n",
    "    dbins = numutils.logbins(1, len(sums), N=n_bins)\n",
    "    spans = list(zip(dbins[:-1], dbins[1:]))\n",
    "    s = np.array([np.nansum(sums[lo:hi]) for lo, hi in spans])\n",
    "    n = np.array([np.nansum(counts[lo:hi]) for lo, hi in spans])\n",
    "    return dbins, s / n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-05-08T02:14:43.343690Z",
     "start_time": "2018-05-08T02:14:39.488036Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf4AAAKZCAYAAABZdklaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XdYVvXjxvH3YYMigrgQHLly5sByj8yy1LQs08qFiqiZK5tqmuYE9545cpSZIzVHZYY4UcuVOULFHAGKgw3n94ff+DUciMADD/frurwuz3nO+XA//uHNeZ7z+RzDNE1EREQkd7CxdAARERHJOip+ERGRXETFLyIikouo+EVERHIRFb+IiEguouIXERHJRVT8IiIiuYiKX0REJBdR8YuIiOQidpYOkBk8PT3NkiVLWjqGiIhIlggNDY0wTbNgWo61yuIvWbIkBw4csHQMERGRLGEYxrm0HquP+kVERHIRFb+IiEguouIXERHJRazyO34REcm+EhMTCQ8PJy4uztJRchwnJye8vb2xt7dP9xgqfhERyVLh4eG4urpSsmRJDMOwdJwcwzRNIiMjCQ8Pp1SpUukeRx/1i4hIloqLi6NAgQIq/YdkGAYFChR45E9KrKr4DcNoZRjG3OjoaEtHERGR+1Dpp09G/LtZVfGbprnBNE1/Nzc3S0cRERHJlqyq+EVERNJiwIABTJ48OXX7ueeeo3v37qnbgwYNYvTo0bzyyisA7Nixg5YtW951rJIlSxIREcH169eZOXNm5gbPACp+ERHJderVq0dISAgAKSkpREREcOzYsdTXQ0JCePrpp1m9enWax1Txi4iIZFN169Zl9+7dABw7dozKlSvj6urKtWvXiI+P58SJE3h4eFC5cuX/nBsZGcmzzz5LpUqV6N69O6ZpAvD+++9z5swZqlWrxuDBg+nUqRNr165NPe+NN95g3bp1fPbZZ7Ru3ZrGjRtTtmxZRowYkXrMsmXLePLJJ6lWrRo9e/YkOTk5w9+7pvOJiIjF9F97lMN/3MjQMat55WNym/8W9t95eXlhZ2fH+fPnCQkJoU6dOly8eJHdu3fj5uZGlSpVcHBwuOu5I0aMoH79+gwbNoyNGzeyYMECAMaOHcvRo0c5fPgwAD/++COTJk2iTZs2REdHExISwuLFi1m2bBn79u3j6NGjuLi4UKtWLVq0aEGePHlYtWoVu3btwt7ent69e/P555/TqVOnDP33UfGLiEiuVLduXUJCQggJCWHgwIFcvHiRkJAQ3NzcqFev3j3P27lzJ2vWrAGgRYsWuLu73/W4Ro0a0bt3b/7880+++uor2rZti53dndpt1qwZBQoUAODll18mODgYOzs7QkNDqVWrFgCxsbEUKlQoI98yoOIXERELetCVeWb663v+I0eOULlyZXx8fAgKCiJfvnx07do1Q35Gp06dWLZsGStXrmTRokWp+/89Lc8wDEzTpHPnzowZMyZDfva96Dt+ERHJlerWrcs333yDh4cHtra2eHh4cP36dXbv3k3dunXveV7Dhg1Zvnw5AJs3b+batWsAuLq6cvPmzX8c26VLl9TZAxUrVkzdv23bNqKiooiNjWXt2rXUq1ePpk2bsnr1aq5evQpAVFQU586l+Wm7aabiFxGRXKlKlSpERERQu3btf+xzc3PD09Pznud9/PHH7Ny5k0qVKrFmzRqKFy8OQIECBahXrx6VK1dm8ODBABQuXJgKFSr85xOEJ598krZt21K1alXatm2Lr68vFStWZNSoUTz77LNUrVqVZs2acenSpQx/38ZfdyNaE19fX/PAgQOWjiEiIndx4sQJKlSoYOkYWSImJoYqVapw8OBB/lpc7rPPPuPAgQNMnz49XWPe7d/PMIxQ0zR903K+rvhFREQywfbt26lQoQJ9+/YlO60oq5v7REREMsEzzzxz1+/ou3TpQpcuXbI+0P/oil9ERCQXUfGLiIjkIip+ERGRXETFLyIikotYVfEbhtHKMIy50dHRlo4iIiLZVGRkJNWqVaNatWoUKVKEYsWKpW4nJCSwdu1aDMPg119/TT0nLCwMwzAYMmRI6r6IiAjs7e156623LPE20s2qit80zQ2mafpnp2kTIiKSvRQoUIDDhw9z+PBhAgICGDBgQOq2g4MDK1asoH79+qxYseIf55UqVYqNGzembn/55ZdUqlQpq+M/Mqsq/sxw8+dN/DHfDzMl4x+NKCIi2cutW7cIDg5mwYIFrFy58h+vubi4UKFCBf5aIG7VqlW0a9fOEjEfiebxP0D8xeNc/2kRGDYU7ToXw0a/K4mIZJSgoCBOnjyZoWOWL1+eQYMGpevcdevW0bx5c8qVK0eBAgUIDQ2lZs2aqa+3b9+elStXUrhwYWxtbfHy8uKPP/7IqOhZQi32AJ4vvINn66Fc37mAK8sHYI1LHIuIyB0rVqygffv2wJ2S//fH/c2bN2fbtm2sXLmS1157zRIRH5mu+NOg4EsjSIm7RdSWSdg45aXQK59aOpKIiFVI75V5ZoiKiuL777/nyJEjGIZBcnIyhmEwYcKE1GMcHByoWbMmQUFBHD9+nPXr11swcfqo+B8gJcUkKiaBwh2CSIm/TcSG0RiOeSjY6kNLRxMRkQy0evVqOnbsyJw5c1L3NWrUiJ9++in1CXxw55eVRo0a4eHhYYmYj0wf9T/AB5tOUGvKT5yKuE3RzjNxq/MGf67+iMitUy0dTUREMtCKFSt46aWX/rGvbdu2//m4v1KlSnTu3Dkro2UoPZb3AQ5cuM7z8/ZiGLC5+1PU8MpL+Ix23Az9mqJ+83Fv1C1Dfo6ISG6Rmx7Lmxn0WN5M5uuTn1196+Fib0vjWSF8d+YaxXqtIE+V5lxa1IPo3SsePIiIiEg2oeJPg3IF8xLStz6lPFx4Yf5eVh+LxKfvV7iUb8jFuR25eXCdpSOKiIikiYo/jbzcnNjZpx5PFXen/bJQZu6/gk//DTiX8iV8RjtuHdlq6YgiIiIPpOJ/CPmd7dnaszatKham79dHGfHjRXwGbsLBqwIXprbh9smfLB1RRETkvlT8D7Bz506GDh1KQkICAM72tnzV2Re/J30Yue0UfTZfwHvQFuwLFOfCxBbEnt1v4cQiIiL3puJ/gPPnz7N582beeustbty4AYCdrQ3z2z3BB03LMHfPeV5fe54iA7di6+rJucDniDv/i4VTi4iI3J2K/wHefPNNRo0axZEjR+jWrVvqmsyGYTD6hQpMbl2JNUcu03J1OB79vsXGwYVzE5oRf/k3CycXEZF7MQzjH6sGBgYGMnz4cABmz57NkiVL7nv+8OHDCQwM/M/+sLAwKleuDMCOHTto2bJl6mtDhgyhefPmxMfHk5CQQP/+/SlTpgxly5aldevWhIeHZ8A7ezAVfxo0b96cGTNmEBERQdeuXTl+/Hjqa/0aPsbnb1Rn1+9RNP3iEnl6bwTT5Ny4piT8GWa50CIick+Ojo6sWbOGiIiI/7wWEBBAp06dMvTnjRo1il27dvH111/j6OjIhx9+yM2bNzl58iSnTp2iTZs2vPzyy1nyPBgVfxrVqFGDhQsX4ujoiL+/Pz/99P838r1ew5sN3Z7kVMRtGq6OwMZ/Aynxtzk3rimJURctmFpERO7Gzs4Of39/Jk2a9J/X/n4137hxY/r160e1atWoXLky+/btSz3u559/pk6dOpQtW5Z58+bd82cFBQWxefNmNmzYgLOzMzExMSxatIhJkyZha2sLQNeuXXF0dOT777/P4Hf6X1qr/yGUKlWKhQsXMnDgQAYNGsTgwYN59dVXAWj+eCG+D6hDi/l7qf/VdbZ0XI3TZ605N/4ZSn64E7t8BS2cXkQk+7n8eX/izh/O0DGdilejyBuTH3hcnz59qFq1Ku++++59j4uJieHw4cPs3LkTPz8/jh49CsAvv/zCnj17uH37NtWrV6dFixb/OXfXrl2cPHmS0NBQ8ubNC8Dp06cpXrw4+fLl+8exvr6+HDt2jKZNm6b1raaLrvgfkqenJ3PmzKFevXqMGzeOKVOmkJKSAsBTJdwJfqsejnY2NNwQT0S7z0mMPMf5Cc+SfPuahZOLiMjf5cuXj06dOjF16v2fvdKhQwcAGjZsyI0bN7h+/ToArVu3xtnZGU9PT5o0afKPTwP+UqZMGUzTZNu2bRn/BtJJV/zp4OzsTGBgIIGBgSxdupTLly8zfPhwHB0debywKyF969N83l6abklhbcv5lFjXhfNBL1B88FZsnV0tHV9EJNtIy5V5Zurfvz81atSga9eu9zzGMIy7bt9r/98VLlyYzz//nKZNm+Lh4UGTJk0oXbo058+f5+bNm7i6/n8nhIaG/uNmwMyiK/50srW15d1336Vfv35s27aNPn36pP4W6J3fmZ196uLr40aLn1w51nQasb/v58LkF0lJiLNwchER+YuHhwft2rVjwYIF9zxm1apVAAQHB+Pm5oabmxsA69atIy4ujsjISHbs2EGtWrXuen65cuVYs2YNb775JocPHyZPnjx07tyZgQMHkpycDMCSJUuIiYnh6aefzuB3+F8q/kdgGAYdO3ZkzJgxHD9+HD8/v9TpGB4uDmzrWZsXHi9Eu4Ne7H1qLDEnfyR85muYyUkWTi4iIn8ZNGjQXe/u/4uTkxPVq1cnICDgH78gVK1alSZNmlC7dm2GDh2Kl5fXPceoVasWixYt4sUXX+TMmTOMGTMGJycnypUrR9myZfnyyy/5+uuv7/qpQUbTY3kzyOHDhxk0aBA2NjZMmjQpdR5nYnIKPb74mcUHwplZdA+NjozGrV4nvLovwrDR710ikvvkpMfyNm7cmMDAQHx90/TE2yyhx/JmE9WqVWPBggW4uLjQs2dPduzYAYC9rQ2L2lfj3Sal6X2pNt+X6kn0riVcWTEwS+ZrioiI/J2KPwOVLFmSRYsWUbZsWQYPHszKlSuBO18JjGtZkcBWFel74wV+KPIaUVunELF+lIUTi4jI/ezYsSNbXe1nBBV/BvPw8GD27Nk0atSIwMBAJk6cmDrdb1Dj0ix+vTr9k1/nJ/fm/LlmGFHbpls4sYhI1tMnnumTEf9uOaL4DcNoYxjGPMMwVhmG8ayl8zyIk5MT48aNo3379ixfvpz33nuPuLg7d/N38vVhXbeneMexF3vz1uXysr5Eh3xu4cQiIlnHycmJyMhIlf9DMk2TyMhInJycHmmcTL+5zzCMhUBL4KppmpX/tr85MAWwBeabpjk2DWO5A4GmaXa733GWuLnvXpYvX86kSZOoUqUKQUFBuLu7A7A7LIqX5gUTFDGMaglHKd5vLa7VMn/+poiIpSUmJhIeHp56QSRp5+TkhLe3N/b29v/Y/zA392VF8TcEbgFL/ip+wzBsgd+AZkA4sB/owJ1fAsb8awg/0zSv/u+8IOBz0zQP3u9nZqfiB/juu+8YNmwYhQoVYurUqfj4+ABw/PJNXpr9HaPDB/O4GU7JwVvI83hDC6cVEZGcJlvd1W+a5k4g6l+7nwROm6Z51jTNBGAl0No0zSOmabb815+rxh3jgM0PKv3sqGnTpsyaNYsbN27QtWtXjhw5AkDFIq5s7/cc40uNJ8z05GxQS2LPHbJwWhERsWaW+o6/GHDhb9vh/9t3L32BZ4BXDMMIuNsBhmH4G4ZxwDCMA3/++WfGJc0gVatWZdGiReTNm5eAgIDU6X4+7s5s6t+C6RUmcSXJiZOjmxF/+TfLhhUREauVI27uM01zqmmaNU3TDDBNc/Y9jplrmqavaZq+BQtmzyfhFS9enIULF/5nul+BPA6sfrsNS6rP4GZ8Ekc+aUJCxHkLpxUREWtkqeK/CPj8bdv7f/us3l/T/Ro2bEhgYCCTJ08mJSWFPI52LOjTlnW1Z5Acc539HzcmITr7fXIhIiI5m6WKfz9Q1jCMUoZhOADtgfUWypLlnJycGD9+PK+++irLli3jo48+Ij4+HntbG8b7v87OBjNxvnWRnUMaE3fzuqXjioiIFcn04jcMYwWwGyhvGEa4YRjdTNNMAt4CtgAngC9M0zyW2Vmyk7s93S86OhobG4N3/Drxy9PTKHTjV7796Blu3Lxt6bgiImIlrOohPYZhtAJalSlTpsepU6csHSfNtm7dyscff4yXlxdTp06lWLE79zluWDKJ0t8NZLd7U1p9soFC+ZwtnFRERLKjbDWdLyuZprnBNE3/v56VnFM8++yzzJgxg6ioKPz8/Dh+/DgArToNILLhh9S59h3zP+5MWFSMhZOKiEhOZ1XFn5PVqFGDhQsX4uDggL+/P8HBwQDU9xtF3JPdaRP1JUEjB3Lk0g0LJxURkZxMxZ+NlCpVikWLFlGyZEkGDhzImjVrMAyD6r1mY1Z+kV5RcxgxYQzBZyMtHVVERHIoFX824+npydy5c6lTpw6jR49mxowZYNhQod8qbB6rx9CoID6YOpf1Ry9bOqqIiORAVlX8hmG0MgxjbnR0tKWjPBIXFxeCgoJ46aWXWLRoEcOHDyfFxo6y72zAqXAZJkV/ynvzv2DhXi3yIyIiD8eqij+n3tx3N3Z2dnz44YcEBASwceNGBgwYQBwOPPbuFlzz5WfB7ZEMWbGdsd+d0qMtRUQkzayq+K2NYRh0796dIUOGsG/fPnr27Em06UyJdzbjYZvIioRRjP1mPwPXHyMlReUvIiIPpuLPAdq0aUNQUBBhYWH4+flxOdkVn35rKRR/ka8IZNaPv9Jx+SESklIsHVVERLI5FX8OUb9+febMmUNsbCzdunXjTKIHxXouo+i1w2x0mcPKg+dpuWAvN+OSLB1VRESyMRV/DlKpUiUWLVpEvnz56NWrFwdjClL49ckUvbCd7UW/5vtTETSZFcLVm/GWjioiItmUVRW/tdzVfz/e3t4sXLiQMmXK8O6777LjZjEKvDCYwkeXsKPCbo5fuUndacGcidD6/iIi8l9WVfzWdFf//bi7uzN79mzq1q3LmDFjWBNZgnx13sBj5xh+qh3GtdhE6k4L5mC4nuwnIiL/ZFXFn5s4OzsTGBhImzZtWLBwEQsvl8KlYlOc1/Xjp2YxONnb0mhmCNt/+9PSUUVEJBtR8edgdnZ2fPTRR/j7+7P+m81MvVAaB69KGMs6sfNFR0q6u/DC/L2sPHTR0lFFRCSbUPHncIZh4O/vz5AhQ9i1/zATzpXByFOAuLlt+eG1ItQp4U6HZQeZsvOspaOKiEg2oOK3En/N9T8WdoUJYWVITk7i2rSWbGr/GC9XKUL/dcd4/5sTWuVPRCSXU/Fbkb/m+p+/ZUNQ2GMkRF3kyrQXWdnucQLqlGDcD6fpuvIwicla6EdEJLeyquLPDdP5HqRSpUosXLiQSEdvpp8vSWzYQS7Neo0ZrR/nk+blWXwgnDaL9nM7Xgv9iIjkRlZV/LllOt+D+Pj4sHDhQm4XrcVnfxTn1i+bufRZT4Y8U5a5r1bl21+v0nT2biJuaaEfEZHcxqqKX/6fh4cHs2fPJqHCi3x9tSjRwZ9x9auh9Khdgq86+/LzHzeoP30XYVExlo4qIiJZSMVvxVxcXAgKCiKlTgA7rnkSueFTIrbPoE2VomzrWZsrtxKoOy2Ywxdz71cjIiK5jYrfytnZ2TFk6FBsXxjB4ZtuXF7al4jdX1D/sQIEv1UPW8OgwYxdbD5xxdJRRUQkC6j4cwHDMOgZ0Js8r8/m91gXLs56nUuh31KpiCt7+zWgrGceWi3cz+yQMEtHFRGRTKbiz0Vav9Ie9+4riEq05+LkFwkL/R4vNyd29qnH848XotdXR3hn/TFSUjTXX0TEWqn4c5n6z7aiYO81JJlwYWILju75gbyOdqztWou36pUk6MezvLrkADEJmu4nImKNrKr4NY8/bSo3eJ7CfdaQ1zaR8MmtCP7uW2xtDKa+VJlJrSvx9dHLPD1rN39qup+IiNWxquLXPP60K1WnJYX9l1PMIYY/Zr7G1199gWEY9G/4GGv+Nt3vnKb7iYhYFasqfnk4XvXbUbjzbCrnvUHEkl7MmT0L0zRTp/tdvZVAvem7OHb5pqWjiohIBlHx53KFmvpT4OWR1MsfRfSGEYwcOZKkpCTqP1aAH3vXJcU0aTB9F7vDoiwdVUREMoCKXyj04ke4P92Llp5XiA2ez8CBA4mJiaGqVz52vVWfAnkceGbOHr799aqlo4qIyCNS8QuGYVCk4zRca7ThzaLhJB3bREBAAFFRUZQq4ELwW/UoXzAPrRbsY8mBC5aOKyIij0DFLwAYNrYU67UclzJ16FPiAjYXD+Ln58eFCxco7OrID73q0vCxAnRecZjhW05imprrLyKSE6n4JZWNgzM+A9bjWOgxBpc+j2vcZfz8/Dh69ChuzvZs7vEUXWr5MGLrb3RacYj4pGRLRxYRkYek4pd/sMtbgOLvfIu9c14+LHueInkMAgICCA4OxsHOhoWvPcHI5uVZFnqRZ+fsISomwdKRRUTkIaj45T8cPEtQfNBmjIRbDC0fzuOlijFo0CDWrl2LYRgMaVaOz9+ozp5z16kzNViP9hURyUGsqvi1cl/GcSr+BN5vf03Sn6f5oPwl6j5Zk1GjRjFnzhxM0+T1Gt5sD7gz17/utGCOXLph6cgiIpIGVlX8WrkvY+Wt1JRiPT4j/tRPDCj3J61atmDevHmMGjWKpKQkGjxWgJ/61MXAoMH0Xfx0NtLSkUVE5AGsqvgl47nVeZ1Cr03g1oHV+Je5Tvdu3Vi3bh2DBg0iNjaWykXzEdK3HkVcHXl2zh7WH71s6cgiInIfKn55oALPD8Lj2f5c2zaVV0vF8sEHH7B792569uxJVFQUJTzuzPWv6pWPlz7bz5zdYZaOLCIi96DilwcyDIPCHYJw9W3LlZXv0KykLRMmTODMmTP4+fkRFhaGZ15HvguoQ/PHCxGw+gjvfXOclBTN9RcRyW5U/JImho0NxXouxfmxJ7k4+3WeLJ6H2bNnExMTQ5cuXQgJCSGvox3rutaiV90SjP/hDO2XhRKbqLn+IiLZiYpf0szGwRmffuuwy1eY85NaUr5oPhYvXoyXlxf9+/dn+fLl2NoYzHi5CoGtKrL6l0s0nbWbP2/FWzq6iIj8j4pfHoqdW2GKD9yImRjHhYktKOTmwvz582nUqBETJ07k008/JSkpiUGNS/Nlp5ocuhhN7anBnLx6y9LRRUQEFb+kg2Oxini/9RXxl08SPv1VnB3sGTduHN26dWPt2rX079+f27dv07aqFz/0rsvN+CTqTA1m5xlN9xMRsTQVv6RL3kpNKdplDrePbePSkt4YhkGvXr0YNmwYBw4coGfPnkRGRlK7hDt73q5PobwONJuzhxUHL1o6uohIrqbil3Rzb+iHZ6sPuf7jfCI3TQDgxRdfJCgoiLCwsNSn+z1WIA8hb9endon8vLH8IDN3hVk2uIhILqbil0dS8OWR5HvqNa5+8R439q8GoH79+syePZtbt27h5+fHsWPH8HBx4Fv/2rSsUJg+a44wattverSviIgFWFXxa63+rGfY2ODV/TOcy9Tl4pyOxJzeA0DlypVZuHAhzs7OBAQEEBISgrO9LV918eXNmsUY+u1J3tlwXOUvIpLFrKr4tVa/Zdg4OOHTby12+b24MPlFEq6eBaBEiRIsXLgQHx8fBgwYwDfffIO9rQ2L21fnrXolmfjjWbqt+pmk5BQLvwMRkdzDqopfLMcuX0GKD9oEKcmcD2xO0s0IADw9PZk7dy41a9Zk+PDhLFq0CMOAqS9VZlizcizaf4HXloYSn6SFfkREsoKKXzKMY9Hy+PRfT2LUBS5MaklKfAwAefPmZcqUKTRv3pwZM2YwcuRIkpKSGNG8PJNaV2LNkcu0mL+P6NhEC78DERHrp+KXDOVSrh7Feq0g9ux+wme2x0xOAsDe3p6RI0fSo0cP1q9fz9tvv82NGzfo3/AxPmtfjR/PRFJ3WjC/R8ZY+B2IiFg3Fb9kuHw121Ck43RuHd7ApcW9U2/gMwyDnj17MmLECA4dOoSfnx+XLl2icy0ftvjX5o8b8Tw19SdCfo+y8DsQEbFeKn7JFB5Ne/1vjv88ItaP+sdrLVq0YObMmURGRuLn58fp06d5uqwne96uj5uTPU/P3s2aXy5ZKLmIiHVT8UumKdh2FG71OvHnmmFc27nwH6/VqFGDefPmAdC9e3cOHTpE+UJ52d23HtWLufHKkgNMD/7dErFFRKyail8yjWEYePnNJ0+V57i0yJ+bP2/6x+tlypRh4cKFFChQgD59+vDtt9/imdeR7wJq06piYfp+fZQPNp7QXH8RkQyk4pdMZdjZ493nS5x8niB8+qvEntn3j9eLFi3KggULqFSpEkOGDGHu3Ll3Fvrp7EvPOiUY+/1pOq84TEKS5vqLiGQEFb9kOltnV4oP3IidW2HOT3yBuPBj/3g9f/78zJgxgxYtWjB37lyGDRuGmZLMrLZVGPV8eZaGhtNywV5uxiVZ6B2IiFgPFb9kCbv8RSgxeBuGnQPnJzQj4cqZf7zu4ODA8OHD6dWrF5s3b2bAgAHExsby0TPlWPRaNb4/HUmTWSFcvRlvoXcgImIdVPySZRwKl6b44G2YSQmcG/8MiVHh/3jdMAy6devGkCFD2LdvH7169eL69et0edKHdV1rcfzKTepN38XZyNsWegciIjmfil+ylJN3JYq/s4XkW5GcG9+M5Fv/nbPfpk0bxo8fz+nTp+nWrRuXLl2iRcXCfN+rLlExCdSdtotD4XoQk4hIeqj4Jcs5l6qJz4BvSPzzLOentCYlIe4/xzRu3Jjp06cTGRlJt27dOHPmDLVLuBP8Vj0cbA0azQzh+1MRFkgvIpKzqfjFIvI83hAv/yXE/hbMH/M6Y6b896796tWrM2/ePFJSUujevTuhoaFUKOzK7rfrU8LdmRfm72Xj8SsWSC8iknOp+MVi3J56jUKvjefGvi+4+uX7dz2mbNmyLFy4EE9PT/r06cPGjRsp5ubMjt51qVzElZc+28/XR7TKn4hIWqn4xaIKPP8O7k17E7lpAlHbZ9z1GC8vLxYsWED16tX5+OOPmTNnDh4u9mwPqIOvd35eXRLKioMXszi5iEjOZFXFbxhGK8Mw5kZH68avnMIwDIq8OZW81Vpxednb3Dy4/q7H5cuXj6lTp9KqVSt2wWecAAAgAElEQVTmzZvHsGHDcLE12eJfm3ol3Xlj+UEW7TufxelFRHIeqyp+0zQ3mKbp7+bmZuko8hAMG1u8e6/AqWRNwme1J/bs/rseZ29vz7Bhw+jduzebN2+mT58+JMfdYnOPp3imrCd+q35m6k9nszi9iEjOYlXFLzmXjWMeig/YgJ1bEc5PbEH8pZN3Pc4wDPz8/Pj00085duwYfn5+RFz+gw3dnuSlKkXot/YYI7ac1Pr+IiL3oOKXbMPOrTDF39kChsG58c1IjLz3R/fPPfccM2fOJDo6mh49enA9MoIvOtakSy0fhm/9jf7rjpGSovIXEfk3Fb9kK45FylJi8FZS4m5wbnwzkm5cveex1apVY+7cucTGxjJ48GCSEhNY0O4J+jcsxdSffsdv1WGSkvVwHxGRv1PxS7bjVPwJig/YSGLUBS5MvvsCP38pXbo0I0eO5Pjx44wYMYLExAQmvliJkc3Ls/hAOK8sPkBcYnIWphcRyd5U/JItuZSrRzH/pcSe2cOlxb3u+519o0aN6Nu3L9u2baNjx46cOnWKIc3KMe2lyqw7doXn5+0lOjYxC9OLiGRfKn7JtvLVaotn62FEB39G1NYp9z22c+fOTJ06lejoaLp27cqWLVt4q34pPn+jOsG/R9FoZgiXbtz7kwMRkdxCxS/ZWsE2H+Naow1XVgziRuja+x5bt25dVq5cSYUKFfjoo49YsGABr9fwZmP3JzkdcZu604L57c9bWZRcRCR7UvFLtmbY2FCs51KcSvlycVYHYn4Lvu/x7u7uzJw5k+eff55Zs2axdOlSni1fiB2963I7IZkGeqyviORyKn7J9myc8lJ84EbsPUtwflIr4sKP3vd4BwcHhg8fTrNmzZgyZQorV67E1yc/O3vXJTHZpOWCfVzXd/4ikkup+CVHsHP1pPg7W7BxcObCpJb3neYHYGtryyeffELjxo0JDAxk2rRplCuYh6+6+HLqz9s8P28vUTEJWZReRCT7UPFLjuHgWQKf/utJunGVC1Pa3HeaH9xZ4nfs2LG0bduWxYsX8/HHH1O/hBtfdKrJwfBoGkzfRfj12CxKLyKSPaj4JUdxLuV7Z5rf6d38Mb8LZsr9F+ixs7Pj/fffp0+fPmzevJl+/frxXGk3vvV/igvX46g3fRcnr+qGPxHJPVT8kuPkq9WWQu3GcWPvKq6sGvzA4w3DoGvXrgwfPpzQ0FDeeustfAs7sqN3HWITk6k/fRcHLlzPguQiIpan4pccqcALg/Fo1peobycSuTkoTee0bNmSsWPHcvz4cQICAngsL+x6qx55HW1pMiuE7b/9mcmpRUQsT8UvOZJhGBR+fRL5ar3KlZXvEL17eZrOa9KkCUFBQYSFheHv74+7Eceut+pT0t2FFvP3seaXS5mcXETEslT8kmMZNrZ4+S/B5fFGXJzXhVvHtqfpvHr16jFlyhQuXbqEv78/NrHX2dmnLjW93Wi3NJQvDv+RyclFRCxHxS85mo2DEz5vr8Wx6OOET32Z+IvH03Ser68v06dPJzIykh49enA76ipb/GtTt6Q77ZeFMvHHM5mcXETEMlT8kuPZ5slP8UGbMBxdOD/5RZJuRabpvCeeeILZs2cTExND9+7duXD2N7b41+blKkUZtP44i/adz+TkIiJZT8UvVsHewxuft9eSFHWB8GltSYmPSdN5FSpUYM6cOdjZ2dGjRw9Cdu5gVceaPF3Gkz5rjhCqu/1FxMqo+MVquJSpjVf3z4g5uZPzk1uREp+2NfnLlCnD0qVLKV++PB988AHfbt7E529Up4CLA/Wn72L5wfBMTi4iknVU/GJV3Op0wKvHEmJO7ODC1JcxU5LTdJ67uzszZszA19eX4cOHE7J9EwcGNOSpEu50WnGYb3+9/xLBIiI5hYpfrE7+em9StPMsbh/dyp9rR6T5PGdnZyZNmkS9evUYPXo0Ids3scHvSaoUceWlRftZeuBCJqYWEckaKn6xSu5N/MnfoCsR60Zy8/DGNJ/n6OjIhAkTqFu3LqNHj+bH7d+ytWdtav/vyv+zfSp/EcnZVPxitYp0nI5TieqEz2hHzJm9aT7P3t6e8ePHU6tWLYYPH86WtV+yxb82z5T1pMeXP7PzTNpmDYiIZEcqfrFaNo4uFB+4CTu3IlwIeiHNc/wBnJycmDJlCs888wyTJ09m25bNrO7sSykPF9ovCyUsKm2zBkREsptsX/yGYVQwDGO2YRirDcPoZek8krPY5S9CicFbwc6ec4HPkRiZ9rn59vb2fPLJJ6k3/K1ZuYwvOtXgRlwSlSfsYNWhi5mYXEQkc2Rq8RuGsdAwjKuGYRz91/7mhmGcNAzjtGEY799vDNM0T5imGQC0A+plZl6xTg6FS1PinS2kxN3k3ITnSL6d9rn5Dg4OTJo0iWeffZYZM2awfMpoQt+uQ/VibnT4/KBu+BORHCezr/g/A5r/fYdhGLbADOB5oCLQwTCMioZhVDEM45t//Sn0v3NeBDYCmzI5r1gpp+JP4NNvHQlXzxA+q32ap/nBnbv9R40aRf/+/fnhhx9YNnMiW/yfoklpT/xW/cz3pyIyMbmISMbK1OI3TXMnEPWv3U8Cp03TPGuaZgKwEmhtmuYR0zRb/uvP1f+Ns940zeeBNzIzr1i3PI83ominGdw+soUrq957qHMNw+DNN98kICCAzZs3M3ViIKverEa5gnlou/gAxy/fzKTUIiIZy84CP7MY8PfPR8OBp+51sGEYjYGXAUfuc8VvGIY/4A9QvHjxjMgpVsi9cQ/iLvxC1LdBOPlUIX/9zg91vp+fH7du3WLp0qXcuHGDDQM/pO70EBrO2MWGbk9Sp6RHJiUXEckYlij+h2Ka5g5gRxqOmwvMBfD19TUzN5XkZEU6TCT+j+NcWuSPQ5HyuJSpneZzDcOgX79+uLm5MX36dOLj49nYawDtv/yNp2ftZtkb1Wlb1SsT04uIPBpL3NV/EfD527b3//aJZAnDzh7vPl9g5+5N+NSXSIx6+LX4O3fuTP/+/dm9ezdzJ4xk11t1qV7MjVeXhLJEN/yJSDZmieLfD5Q1DKOUYRgOQHtgvQVySC5ml7cAPv3XkRJ/iwuTXyQl7tZDnf/Xd/6DBg0iNDSUDV98zpYetWhS2pOA1b9w9NKNTEouIvJoMns63wpgN1DeMIxwwzC6maaZBLwFbAFOAF+YpnksM3OI3I2Td2WK9VpJ3PmfCZ/1+kPd6f+Xl156iYYNGzJr1iw+GTaEpR2eIJ+TPe2WhnIzLikTUouIPBrDNK3n63DDMFoBrcqUKdPj1KlTlo4jOUTUtulcXtYXj2f7U+SNSQ99vmmaLF26lKlTp9KxY0cqPf86z83dQ4XCrgS2qkjzxwtlQmoRkf9nGEaoaZq+aTk226/c9zBM09xgmqa/m5ubpaNIDuLR7C08nu1H1NbJRG2b/tDnG4ZBx44dadeuHUuXLuXazzv4tkdtbsYn8fy8vYT8/u8ZrSIilmNVxS+SXoU7BJG3Wisuf97voZ7m9xfDMBg4cCANGjRg/PjxxJ/ay/HBjSmY14F3vznO1ZvxmZBaROThqfhFAMPGFu9ey3EqXo2Ls9oTd+7wQ49hZ2fHmDFjUtf2/2btV4x5oQK7z12j7rRgYhL0nb+IWJ6KX+R/bJzy4jNgAzYu+Tk/qSWJUQ8/y9TJyYmpU6fSsGFDAgMDSTy0ic3dn+RMZAzPz9uru/1FxOKsqvgNw2hlGMbc6OhoS0eRHMre3YviA74hJTaaC5NaPvQ0P7jzVL/Ro0fTokUL5s+fz+5Vswls+TjHr9zimTl7CL8emwnJRUTSxqqKXzf3SUZwKv4ExXp/QdyFXwgb+zRJ1y8//BhOTnz88cf07NmTDRs2ELF1Id/3fJIbcYm8teYI1jSbRkRyFqsqfpGM4vrE8/i8/TXxF48RNu5pkmMe/lMkwzDo0aMHffv2Zdu2baxbNIMRz5Vn3bErvP31UZW/iFiEil/kHlxrvEjxgRtJuHKKi7M6pGuBH7izvG/nzp1Zu3Yt3pf3M7BhSabvCqP3V0dISVH5i0jWUvGL3EeeCo0p2nEGt37ZzJWVg9M9Ts+ePalduzZBQUHc2jSdwQ1LMnv3OT7c9GsGphUReTAVv8gDuDfxx6PZ20RtmcS1HxekawwHBwemTZvGO++8Q2hoKA1tz9GzTgnG/XCa709FZHBiEZF7s6ri1139klkKdwgiT5XnuLS4F7d//TFdYxiGwWuvvUa5cuUYN24crfNdoXQBF56ft5fAH87oO38RyRJWVfy6q18yi2Frh3evlTgUeozwaW2Jv3QyfeMYBkFBQZQqVYqxn45ioq8tLSsWYvA3x1m8/+EfDywi8rCsqvhFMpNtnvwUH/ANGDaEfdqA2N9D0zVO0aJFGTlyJHny5OGTD9/hJfvfaPxYfvxX/8zmE1cyOLWIyD+p+EUegkPhMpQcEoyNowvnxjYm5rfgdI1TokQJvv76a3x9fZk8aRItEw9RqbArr39+iE0qfxHJRCp+kYfkWKQcJYeEYOdWhPBZHUi+fT194zg6MnPmTJo0acKaL1cxrKpB4bwOtJi/jw83nSBZU/1EJBOo+EXSwd7di2I9Pyfp+iUuL++f7nFsbGx45513KFasGJNHDWX1q2VxtrdhzHenWXHo4Z8VICLyIFZV/LqrX7KSc+kn8Wz5AdHBi7lxYE26xylcuDATJ04kKSmJTz4YRGjn0lQvlo8B644RFhWTgYlFRKys+HVXv2S1gq2H4lSqFhfndkr3ND8Ab29vpkyZwvXr1+nduxdj67qRlGLSYv5eLt+Iy8DEIpLbWVXxi2Q1w84Bn35rsXP15NyYxkRsGJPusXx9fZk7dy4JCQkM7deTCU86EHYtlleXhBKflL7lgkVE/k3FL/KI7N29KP3pUVx923L162HEhh1M91glS5Zk1apVFC1alDWzA+lX053g36N4IvBHzkTczsDUIpJbqfhFMoCNU16Kdp2DnWshzgc2J/6P9K/B7+3tzdixY7lx4wb7Zg9l+FNuXLwRR51pwUTeTsjA1CKSG6n4RTKIXd4ClPjgB8DkjwV+mCkp6R7r8ccfZ9myZdja2rJn0Ri+aluKqJhE3vvmRMYFFpFcScUvkoEci5SjcPsgYk/v5vrO9D3Q5y8lS5ZkxowZxMfHs3jiKAY0KMmCfecJ2nEmg9KKSG6k4hfJYG71OuLyeCOufPEe0btXPNLDd8qUKcPgwYM5efIkPmHf83KVIry/8QSzQ8IyLrCI5CoqfpEMZhgGRbvMwS5vAS7Ofp3ITeMfabzmzZvz0ksvsWTJEoqELqVpKVd6fXWE9UcvZ1BiEclNrKr4tYCPZBeORctTeuxJ8tXuwNUv3id676p0j2UYBh9++CH9+vXjwL69eATPpaKHHe2XhfLZvgsZmFpEcgOrKn4t4CPZiWFjg1f3RbiUa8Af8zoTcyok/WMZBh07dmTcuHGcPn2Kp//cRgVPJ3p8+TMHw9P3rAARyZ2sqvhFshsbe0e8+32NvYcPF6a0IeHq2Ucar3HjxnTu3JndwTvx2TsXB5JoMCOEHacjMiixiFg7Fb9IJrPLWwCfgRsxU5I4P6llup/m95devXoxaNAgwn8/w9iykZR0d+aF+Xv57c9bGZRYRKyZil8kCzgWKYfP21+TcOU04TNexUxJ/xK8NjY2dOjQgfr167N0/hzeLxGFYRj4f/kLN+OSMjC1iFgjFb9IFsnzeCOKdpnN7WPbufx5f1ISHu3hO2PGjKFWrVpMmTCGIVUMgn+P4rm5e7gem5hBiUXEGqn4RbKQe0M/3Jv25tr26Vxe0vuRxnJ2diYoKIiCBQvyy4YlDGvgxe5z13h6VggxCbryF5G7U/GLZLGinWbg0Xwg139aROTmoEcay8nJiaFDhxIWFsauGR/QzzuSQxdvkOeDzVyL0br+IvJfKn4RCyj44hCcy9bjysp3uH3ih0caq27dukycOBGA4M+nU8fxTwD6rzv2yDlFxPpYVfFrAR/JKWzzuFPi3e3YuRfj0qKeJEY+2kI8derUYcOGDbi6upKwdSa9q7iwNDSc0dtPkZic/ocFiYj1sari1wI+kpPYODjh3XsVSTeu8PuoesRfOvlI4zk4ONC6dWsAzqwKpH7+OD7a/Csjt/2WEXFFxEpYVfGL5DQu5epR4oMdmEnxnBvbhMSo8Ecar3fv3nz66adERUbiuvszXqpYkJHbTjFvz7kMSiwiOZ2KX8TCnEtUp8R735ESd4vzk1qSEpf+hXgcHBx47rnnCAwM5MqVKzS1CwPA/8tf2Hf+WgYlFpGcTMUvkg04eVfGu88XxF84wvkprYkNO/hI49WvX5+KFSuyaFogI8vdueflqSnBHArX/S8iuZ2KXySbyFu1OV7dFhD7WzC/D69F7LlD6R7LxsaGGTNmUK5cObZ+PpcAt7MYKUn0/foIpmlmYGoRyWlU/CLZSP4GXSg99iSGvRNXv/zwkZb2dXV1Zfz48ZQrV47965fSLe9pdoVd483PD5GkO/1Fci0Vv0g241CwJIXbjeP2kW+J/HbSI43l7e3N/PnzqV69Ohf2f8f7tT1Zfugiz8/bqyt/kVxKxS+SDXk0e4u8T7zA1S/e5drOhY88XkBAADdv3mTbuL7U5ALbT0UwcP0xlb9ILqTiF8mmivkvxalULa6uHEz8H78+0lg1a9Zk7ty55MmTB9vgxbRzu8Tknb/zxeE/MiitiOQUKn6RbMo2rwdFO83ETEkibHQDYk7vfqTxKlWqxJQpU0hOTib+4CYe93Sm37pjHLt8M4MSi0hOoOIXycacS9WkxHvfkRxznXNjmxB34cgjjVetWjUmTZpEeHg4lX9dRdSfV6k99SeOXrqRQYlFJLtT8Ytkc86lfCk98mcAfh9Zh/jLpx5pvAYNGtChQwfO//oLz4V/hb2ZRMsF+zh/LSYj4opINmdVxa+H9Ii1cixWkVJD92DY2nPmvXJc2zHvkcYbNGgQU6dO5dLFC9T+bQXhV6N4IminpvmJ5AJWVfx6SI9YM6cS1fDyu1P4lxb5P/Ld/k899RT9+vXj6rlTVAsZT/TtWKoG/ajyF7FyVlX8ItYuX61XKDftKg5FynF19YeYyUmPNF7Hjh2pWLEiAIXC93Diyi2enbuH5BRN8xOxVip+kRzGLl9BCrUbR3L0FcJGNyQl/tG+m587dy5Vq1bF+/fvKBF1hB9OR+L8/kZd+YtYKRW/SA7kWr0VDkUfJ/b0bn71z/NIj/N1cnJi4sSJ1KhRA88jayh2ZhuJySYvLtyvBX5ErJCKXyQHMmxseeyTQ7iUbwjA+aDnSbp+Od3j5c+fn2nTplGuXDm8r+wnf+wVNv96laAdZzMqsohkEyp+kRzKxsGJEh/soNBrE4gPP0r47NcfaTxHR0cmTpxIfjc3qhxbim1iDIO/OU7I71EZlFhEAE5cucmZiNsW+/kqfpEczDAMCjw/iHy1XiXmxA9cmNKGxOuX0j1ekSJFmDp1KvFxcTQ+sxyH2CjqTd/F7fhHu4lQRP5fxfE7KDPme4v9fBW/SA5nGAZePT7Do/lAbh35lj/md32ku/3LlSvHmDFjuH7lItWOLAIzhbwfbrboFYqIZBwVv4gVsHF0oUiHIAq0eJ/bR7YQPv0Vkm5GpHu8p59+mjp16pAce4uaO0dipCRRcfwOXfmLWAEVv4gVKdh6GAVfHsnNg+u4MKnlI40VGBiY+vcaP31KYmICeT/crDn+Ijmcil/Eihg2Nni++BEezd4m9sxebuz7Mt1jOTo6sm3bttTtKnsmg5mC36rD3IzTlb9ITqXiF7EyhmFQ6NUxOJWsSfiMdlyc14WUhNh0jeXu7s7WrVsBsE+8TbXgsSwPOUnH5QczMrKIZCEVv4gVsnF0oeRHwXi2+pDo4MX81rcwt0/+lK6xPDw8WL9+PQC2KYmUP7yIdceuMHf3uYyMLCJZRMUvYqVsHJwo9MqnuDftQ0rcTSI3TUj3WF5eXuzbtw8Ap9hI8v95jJ6rf+GNZQe1up9IDqPiF7FyRTtNJ39jf279sokrX36ImZK+NfhtbP7/v4vSx1dT5NxOVoSeY9JOre4nkpOo+EVygUKvfEreqi8Q+c0YLkx+ETMlOV3jLF++PPXvxcJ+oODFfQxaf5yIW/EZFVVEMpmKXyQXsHP1xOftNeSt/iK3ft7Ib/2LkRh18aHHKVeuHPv27WPixIkA+JzdhvvVIxT8eCuXbsRldGwRyQRWVfyGYbQyDGNudHS0paOIZDuGrR3eve5csSdHX+H3EU+SHHvzocexsbGhYcOGTJo0CYDHTqzBSElm9PZTmuMvkgNYVfGbprnBNE1/Nzc3S0cRyZZsHPNQ4r3vAEi6/gcnA/Kle6pfgwYNGDp0KACN4w8zfVcY9u9+w5Wb+thfJDuzquIXkQfLU/Fpig/emrodsX5Uusdq3bo1NWrU4Ma+b/A5tQnThCLDt+o7f5FsTMUvkgvlrdyMx+fF4urblogNo7kwpQ0piekr69GjRwNQ6I/9PLFrPEZKEgU/3kpCUvpmD4hI5lLxi+RSNg5OFOk4DYCbB9fxa3cnzKSEhx7H09OT77+/84hRu6RYngi5s17A8/P2ao6/SDak4hfJxezzF6Xc1Mup2ye6ORL7+4GHHidfvnz88MMPANgmJ1A07Ee+P3WVoiO2PeBMEclqKn6RXM7OrTClx53EoWh5AH4fXitdy/u6urpSunRpALzO7aDS/hlcuRnPtJ9+z9C8IvJoVPwigmORcpQcsosCLd4H4Nzohtz8efNDj7Nq1SpeaHHnccBOsVGUObKct9ce5YvDf2RoXhFJPxW/iABgl7cAhV4dTekxJ7B1Lcilhd3T9R39sKFDGD9+PABuUafwuPILry3ex4krD79mgIhkPBW/iKQyDANHr8fxfPEjkq7/wfnxzUi89nBX63Z2djRp0oSKFSsCUOrXr6m6ZyIVx+8gRQv8iFicil9E/iN/vU44lajO7ePfcap/MY53Nog7/0uazzcMgyVLlrBy5UoA7BNjyB9xAtvB33AjLjGzYotIGqj4ReQ/bPO489gnBynWa0XqvrNDnyDx+qWHGqdMmTKpfy997Au8T3+L+3vriLz98NMGRSRjqPhF5J7carenzPhTqdvnxzcjKfrKQ42xZ88eOnToAEDhi3spf3gRnsO2cCs+KUOziuQUHld+Ie/1MIv9fBW/iNyXQ+EyVFxsUqznMuIvHuPc+GdIibuV5vPt7OwYOHAgDRs2BMDl9hVcboTj+uFmPg8Nz6zYItlWqV+/pvzPiy3281X8IpImbnXfwKV8Q+LDj/JrT1eivpuZ5nMNw2DixIl0794dgAqHFuAYG8Wbyw/xR7Qe5yuSlVT8IpJmxQI+x9W3LQCXl/ThuJ8DKfG303z+G2+8kfr3yvumYR9/g2KfbCMpWev6i2QVFb+IpJm9hzc+fVdTrOeyOzuSE/nVPy8xv+1K0/murq7s3bs3dbvqnkkUO7udOtOCMyOuiNyFil9EHlq+Oq9TJvAsGHf+Cwn7tD63jmxJ07m2trZs3vz/qwIWubCL0LAIeq3+hbjE5EzJKyL/T8UvIg/NMAwcCpaiwvw4DAcXAM4HNufinE7EnTv8wPMLFixIcHAw1apVA6BG8Gg+27YPl8Hr9EQ/kUym4heRdDPs7Hl87i0KvDAYgOiQpZwdVj1Nd/07OTkxf/58ChYsCECl0NnUCB7NhmNXiE/Slb9IZlHxi8gjMQyDwq+Nx6l4tdR9v/b2IDkmOk3njxgx4h/br03fiNN7m3TDn0gmUfGLSIZ4bOQh3Jv2ubORnMjJXvmJ3vvFA8+rVavWP+72r3RgFoUvhGD/7ka2//ZnZsUVybVU/CKSYYp2mk752f9/pX9x5mvc/HkTiVEX73mOYRgMGDCAmTP/f10A77PbsE2ModmcPSTqyl8kQ6n4RSRD2Trno0xQGLaud767vzCxBacGeD/wvFq1atG/f//U7WohE3CL+BWHwd/oqX4iGUjFLyIZzsGzBOWnX8W771ep+8Knt8NMuffVu2EYvPnmm7z22mup+8ocW0XV3UEUHr41U/OK5CYqfhHJNPl8XyZ/ozvL9N7Y/yWnBpUk7vzP9z1n8ODBBAQEpG7bJ94m4nYCZcd8ryt/kQyg4heRTOXlN4/HRv0CQFLUBc4OrUbEhjH3nfLXuXNn3nnnndTtsr8s5Xz4Rbp/cf9fGkTkwXJE8RuGkccwjAOGYbS0dBYReXhOPlXw6r4odfvq6g859W4ZzKSEux5vb29P+/bt2b59OwD5rp2lyt4pLAk5ScXxP2RJZhFrlanFbxjGQsMwrhqGcfRf+5sbhnHSMIzThmG8n4ah3gMePC9IRLKt/A26UHGxeWepXyA5+gonujlyvLOBmZR493Py52f+/Pmp29VCJuDyxWDavz8hSzKLWKPMvuL/DGj+9x2GYdgCM4DngYpAB8MwKhqGUcUwjG/+9aeQYRjNgOPA1UzOKiJZwKFgKcoEhf1j3x/zu95zyl+1atXYvXv3P/ad3r4Ku7dXZVZEkf9r787j66jr/Y+/vjlLTvakSbqvNF1oCy1tKYisKsrSijugCEiheEFQRAVFQUThiuCVzauoKKiXRX4gFL0CcmWTrS2UpXub7kuSNvt6tu/vj6R7lpNkTuZkzvv5ePTxyJn5zsy7nSafzMx3vl9PS2rht9a+DFQfsngesN5aW26tDQOPAOdYa9+31s4/5E8lcCpwPPBF4DJjzKB4PCEiXQuWjGPiT9fu+1z3+p9Zd81oal76XaftA4EAl1xyyUHLZr32M8w3n0pqThEvcqOIjgK2HvB5W8eyTllrb7DWfgP4H+A31tpO3wcyxizq6AewtKpKo32JpLrM4ZOY9qAFX2Dfsp0PXEr9W38hUrPjsPZXXHEFr776KkNG7B8TYM7Lt5B76S8PaysiXRs0V8/W2j9Ya5/pZpsIm3sAACAASURBVP391tq51tq5eyf9EJHUN/GW5WQfedq+z9vu+wLrvjGKSPW2w9qGQiGeW/zXg5ZNXf4A5trFSc8p4pTj8qspy+p5IqtkcaPwbwfGHPB5dMcyEUlDmaOmMf76/2PExb8+aPm6a8bQuOKftG5bcdg2jzzyyEGfx615mq8+/h5Pf7ArqVlFnHDlmI3ceMQa147vRuFfAkwyxkwwxgSB84CnXcghIimk6LRFTP1tKyMW7n/Ov+X20ym/YQZVT91CrKVh3/KysjJ++MMf7vtcsusdlv7nJXzqd68T0yA/It1K9ut8DwOvA1OMMduMMQuttVHga8CzwCrgMWvt4b/S9+14C4wx99fVJTYdqIikloxAJoUnXkzBiRcftLzqiRtZ89V8al78zb5l8+fP58033zyo3exXfsKUm3XbX6Q7xlrv/XY8d+5cu3TpUrdjiEg/hCs2sP3XX6Jlw5uHrRt/4xtkTTgWk5HBnj17+MQnPnHQ+q/e+RCXnjJtoKKKJKy1tZXyy7MA2ju3OsQYs8xaOzeRtoOmc5+IpJfgsIlMuPENjvx9jLzZ5xy0btOPjmfD9VOIVG+juLiYJ58+uN/vr669kPEX/5i3t9UOZGSRHkWjUbcjqPCLSGozGRmMvvpJJv704M5Q4Yr1rLumvZ/wmJHDufHGGw9aX/LBX/n4dfdQ29L5qIAi6UqFX0RSnjGGzOGTGbnoj4etW3mRAeCTn/wkb7311kHrxq95irKFt/H6pkPHERNJXyr8IjJoFH74AqY9aJl8186Dlq+8yLD7bz/FxKO8+uqrB60bv/Zpzv3691lb2YCIeKzwq1e/SHrwFw5n6q8PLuSVj13PqoWZ+KNNLF26lJtvvnnfuqHb3+KLZ53GP156baCjiqQcTxV+a+1ia+2igoICt6OISJJlhHIZf+Mbhy1fe2UJNf+6n7PPPvuwdd+/9mouvf0PVDS0DUREkZTkqcIvIukle+Jxh830B7DzD5ez8iLDM985nV/d8s2D1i1/7F4+cvanaGjQrX9JTyr8IjKoBUvGMe1By8RbDx8HrPpv/0n2I1/i2+ccA+x/ZzqruYrj5p/H/66qGMCkIqlBhV9EPCFz1LT2jn93Hz5e/1Hrf8tD098mz7f/1b7spgp+8OWz+dQ1P+ZPyw6fEEjEqzxV+NW5T0T8BcM48vexTtfdN/U97p783kHLtr3yV2795uXEYp3O+C3iOZ4q/OrcJyLQPujPtActE25667B1hYEID01fxtBgK8fm1wDtV//HHTePnz/2HGsq3ZsuVbxvZYX7fUs8VfhFRA6UdcSxTP1NCyMve/CwdXdMWsFVY8p5aPoyzireBVj+5/bv8dnPfoalazcPfFhJC7sfv8HtCCr8IuJtGcEQhSdeyLQHLcPOu6PTNucN3869U95jYlYjoZZqzr36B1x2xVVUVVUNcFrxuvErf+92BBV+EUkfxWdey9Av/LTTdfn+KDcdsYaHpi/jiNrlvPPW65x55pk0NurWv3iL3+0AIiIDqeTs75A5ajo23My2+77QaZu7prxPbcTP45WjWHzhKIafehElp13G1KlTCQQCA5xYxFkq/CKSdvJmtY/qN22epXnDmxgMG3903EFtCgNRLh3V8ax/2T1899EXOGrKRD4/s4iRn/gPSibPw2TopqkMPp4q/MaYBcCCsrIyt6OIyCCRPbG94BeffR17/tb5YwCA28pWQmwlvA1Vbz9EFZA18XjGf+8ljD84QGlF+s9Tv67qdT4R6auhn7+NI3/Xxrjvvki8eGJC27RseIOW8iVsvOXDVD9/Ly3lS5KcUqT/PHXFLyLSV8YY8AfJmXoKM36+njVfKyXWsLvH7Tb95EQAWta3z/yXO2sBRacsJFAygdDYo5OaWaQvPHXFLyLilMm/2MHke6uY9qDlkiE/Y0NzNr/fMbbH7RqXL2brXZ+i/Aczad32AbHG6gFIK5I4XfGLiHTC+AP480oAuGnRlznvdh+Zvhr+teLvlGU1cuMRa3rcR/kNRwEw5htP48sfuq8/gYibVPhFRHpw5pHDePg7X6Q0J5M3Nn6ZO+64gwtX5BI0cb4/YTXjs1q63X7rLz657+upv24gI5Sb7MgiXdKtfhGRBJx15DCOHVvIVadMYsNf/5tjL7uZsM3g6tpPceGK2TxZOYLqSIB7t07odj+rL89j5UWGppX/h41FByi9yH7GWttzq0Fm7ty5dunSpW7HEBGP21zdzPifvICJRZj96q37lmdmxPjNkcsT3s+kn28hUDwmGRElxay8yOz7etqDztVfY8wya+3cRNp66opf0/KKyEAaNySbd689hVXf+zhfvu5WmnOHA9AW93HpyllctGL2vrb/s2t0l/tZ982xrP36KMIV64m11BMPtxBrqkl6fklPuuIXEXFIPG4Z89lrGLH11U7XlwbauHPyBwnvr+ScGyk86WKCpd0/PpDBQ1f8IiIekpFhWPbALZz4+UvZOPVTh62vimRy4QF3AXqy+6kfsf5bR7DyIrPvT6R2J5tv/zjxVk0eJH2jwi8i4qDhRXn84rqv8ufvL6KuqLMRAA2XrpzV5/2v+/pImlY8z+rL89j+qwv6HlTSlm71i4gk0a3/XMcN/7sagDkv3dxpm2vGrueYvP71TRq58AFyjz4Tf+Hwfu1Hkku3+kVEPO6T04cBMGtkPm+fdAOrZ11yWJv/2jKR9xvzeaJyRJ+Ps+N3l7D26yOoffUh6pc8DkDjB89T+eQP+7xP8SYN4CMikkQzRuRj71wAwJ6mMCU3PkvlyGMZuuPACX0MP9s8CYCnq0aQ5YtxdvEu5pdW8HTVcD5ozOfikVsYmdna4/F2/Oaiw5bZtiYaP3iOrAlzGfKxqwiN6/ujBhn8dKtfRGQALdlSyzefXsGr5XuY+e+f4o+1Jbxtvi/CvVPfcyTH1PubyMjMdmRfkrhUuNWvK34RkQF07NhCXvnah1lT2cjUn+4vAoG2esau+zuFe7qeA6A+FmDhymP43bR3+p1j9aIcAPKPO4/RVzzc7/3J4KFn/CIiLpgyNBd75wI+NK4IgEhmPhtmnMeKOV/tdruIzeC328dxz9Yj2NyS1e8c9W8+wq6Hv9Xv/cjg0atb/caYImAk0AJsstbGkxWsL4wxC4AFZWVll61bt87tOCIiPWpqi1LTEmF0YRa//PcmrnzifUp2LCMjHqE5dzhT3n2w2+0DJo4BMozl+xPWMDbU/YRBPRl+wd3UL3mcMdc8gy8rr1/7ksOlwq3+Hgu/MaYAuBI4HwgCVUAIGAa8AfzSWvuvfiV2mJ7xi8hgFY7GmfLT/2NT9f4Cvvc1wM2Tzmbcur8ltJ/sjCifGrqTM4or+5xl6Lm3U3LWt/u8vRxusBT+54GHgMXW2tpD1s0Bvgy8b639XR/zOk6FX0QGs3A0zrNrKvncg8sIx+Jg42TEo8R9QbCWIRXvMWHNX3u93ynZDdwwYW2vthl95WNkjp1F5vBJvT6eHG5QFP7BSIVfRLzo2qdX8POXygHIiLaRW7eFSR/8T6/2MTm7gUtHbmZ4ZuJvEwCM/8HrROt2kXvUGWQEQ73aVvYbVIXfGPNhYLm1tskYcwEwG7jLWru571GTQ4VfRLzMWssnH1jCMyt2MW7t0zQUTmDC6id7vZ+TCndz2aje/QjPyC5k6n/XEG9txFqrfgC9lAqFvzev8/03MNMYMxO4Fvgt7Y8ATul9RBER6StjDIsXzuPlDXs45ZfthaR66FFg2r/uamjgQ71SW8JrdUO4anQ5s/MTGzI43lx7UPGa+psWqv95DwXHnUegeEwv/ybiht68zhe17bcHzgHutdbeB+hXPRERl5w8sZjtN55OQci/r+gDrJvxRQDaQkU97iNmM/jF1jK+u34am1qy+FdD74r36suyqHz0O6z75li8+OjYi3pzxd9gjPkucAFwsjEmAwgkJ5aIiCRiZEGI2p+cSW1LhIff2c6/1u/mL+/CslNu2tcm2NreL/uoN+/qcj/b27K4sXwaAP8XyuVDBdWcVVLRqyy1L/6GotMW9eFvIQOpN8/4hwNfBJZYa18xxowFTrXWPpTMgH2hZ/wiku6qGtsYetNzhy0fu/YZSncuS2APlvOHbeO4ghqGBCIJH9fJ59ZelArP+Hu81W9M+/0ja+0ua+3PrbWvdHzesrfo720jIiKpoTQ3k79+5VgAirL235zdMnk+y06+MYE9GB6uGMM31h7NbZsSf5Vv5UWGrfd8llhz/6YZluRJ5D3+F4H/Bzxlrd1ywPIgcCJwEfAva+0fkhezd3TFLyJyMGstT76/i//4f+9R2Rjet7x0+xLGrv97Invgm2PXMyuvPuFj6ur/cIPiih84A4gBDxtjdhpjVhpjyoF1tI/m94tUKvoiInI4YwyfOXoEFTd/AnvnAnKCPgCqRh2b6B54cOe4fZ8e2TWqxy3Kbz6OWHMdbbs0hHoq6e1Y/QGgBGg5dBS/VKIrfhGRnr2+qZoT7vk3M968i8zWxH6klwTa2BMJYmm/cp1fspMvDNvR43ajr36S/DmfAiBaX4kvtwSTkX7zxKXCFX8it/pDwFeBMuA94AFrbbTfKZNAk/SIiPTO31dVMDxkOOH2v5Fbv436ookU7l6d4O3/dpeM3MSpRXsSajvk41+n+rm7KFnwPYZ+7id9jT1opULhT+TXrQeBucD7wFnAnf3IllTW2sXW2kUFBQVuRxERGRTOOnIYsycM5flvf5I9w2cRycyjatSxrJy9iOrSaQnt44Ed43pu1KH6ufZXCuuX9X6kQXFGIoV/mrX2Amvtr4HPASclOZOIiAywk44o5g/nzWLZNe0/4lvyRrB58icT3NrwnXXTe3W88I5VbL37M6y/bgrx1sZeppX+SKTw73uBM1Vv8YuISP9ddOwYZo8u5MmL2+8Yx/2ZvHfcNxLadlc4xH1bJ/TqeA3LniS8ay2rL0+PQWBTZWTDRAr/TGNMfcefBuDovV8bYxJ/r0NERAaFTx01Yl/P/0iogPfnXcXOsSdRVzSx2+3erB/ChSvmcPmqWb0+ZrS+qq9xB42mpia3IwAJFH5rrc9am9/xJ89a6z/g6/yBCCkiIgOv5sdncPrkEsJZQ9gx4SOsP/oClp18I4353Y/n3xL3ceGKOb061tqrhrLyIkM83NKfyJKA9HuXQkREEhLwZfDc5R/i4Qtm719oDGuOuYRlp9y0bzKgrly4Yg6/3T6O5Q2JXyOuviybWFNNXyNLAlT4RUSkW587egQ/OH0StT8+46Dl9cU9D+X7cm0JP9+S+JC/AGuuGNKr9tI7vRrAZ7DQAD4iIslhreXfG6tpjsT4xP1vYmIRSncuY8yGZ7vdrsAf4dxh2zixsLpXx5t8TyVkZODPLe5P7JTQ2NjIliv3d2RM5ff4RUREgPahf088opiPTxlK7GfzufTDE6kcfTwbpn2h2+3qogHu3z6BH5VPoTmWeOlZe9VQ1l5Z0t/YKaHx9T+5HQFQ4RcRkT7KyDDc//mZXHPyEdSWHpnQNutbcrli9Sye31Oa5HSpp/pP/+F2BECFX0RE+unn57QP3vP+vKuJG1+P7eMY/rhrLBeumMPVa45K6BgrLzLEGnv3mCBVtGx6m1hL6rz9rsIvIiL9ds3JRxDOKuKdk7/P2yfdkPB2tdEg16ydkVDbNVcWE6neNqje+bfRCBtvmsPWXyQ6CmLyqfCLiEi//fyc6cyfNgwAm+Hv1bZ7IpkJt113zRjWXjWUcNWmXh3DLdbGAWhZ/7rLSfZT4RcREUcsXjiPGcPbe60v//B11BWV0VAwPqFtezvgz/pvTSBau6u3EQUVfhERcdD73z6VptvOJOYPsf7oL7F21kUsO/nGhLa9aMXsnhsdoHndv/sS0RXxFHp1XoVfREQclR30c+tZU/cvMIaVsxf1uJ3F8O1ezPK37d7PUfXULX2JOPBU+EVExMuuO62Mq07cP1tfS96IhLarCIf45toZfD3B3v5VT9xI43v/6FPGgWCjYQBMPHUmt/VU4TfGLDDG3F9XV+d2FBGRtJaRYbj70zOI3zF/37L3jr8moW13RzKpiQYTPtaWO888qKe/jcepfv7elJjwp+Htv7od4TCeKvzW2sXW2kUFBQVuRxEREdpH+tt4w0cBiGTmU37k5xLe9sIVcxJ+7r/2qqFs+a8F1P77j9S/9Ri7/nQVVU/e1KfMjkqhW/x7earwi4hI6hk/JBt75wJOGF9EzdDpLDvlJlbPuoSKUcf1uG1vnvs3Ln+GHfdfuG92v1SY5W/lyhVuRziMCr+IiAyI5y8/nrevOZkL5oyiqWAM28rOoDWr55n4KsIhfrVtfMLH2fXQFf1I2T/h3ZsP+vzoo4+5lKRrKvwiIjIgsoN+jhldwB+/uP/2/Ypjv0b5kZ/tcdvX6voyO5/pwzZ917jiBdZfO5661x/etyz1bvSr8IuIiAteuuKE9i+MoWboDN790Ld63ObHGyfzo/IpCR+j9pUH+hqvT9q2vgdAS/lbA3rc3lLhFxGRAXfyxGJiP9vf4z8azOlxoJ+1zXmsb8lNeGIf4jHq3/oLNh7HRiP9iespKvwiIuKKjAyDvXMBM0fmty8whj3DZva4XW00yKu1PfcNANh23xfYdNsprFoYJB5u7U/cBLTf2H9tU2rPIqjCLyIirpq1t/ADkUB2Qtu8VVeU8P5b1r4KwIbvTu2hpTPe3JLaY8mo8IuIiKvu+tQM7lgwjcmlOewaexJ7hh3NOx++vttt3m/s/XgtkUN63KcrFX4REXFVQVaAa0+dyPThecQCWWya+mni/kx2jflwl9vE+thjv2rxrbSULyHauKevcbuWgoP1dEaFX0REUsKD5x3DMwvn7fu8/YiPsXPsSV22/+baGexqy2RjS2KPBwCqHr+BjTfPY+2VJf3Kuv66Kez53zs7XXda+M2Unj9AhV9ERFJCXsjP2dOG8dyi4/ct2zHhI9SUdP5sfnckk++sn8FN5cl5dr/xJyex/VcXHLa8ZeMywrvWUvFI+yuI9UufYNWlIeLhZgDGxnex5c4zk5LJCSr8IiKSUo4fV0RhVmDf561lPRVRwzsNvX/mv+UX5wBQ+cRNbL7j8GO0rH2Vutf/fNjyjT+ce9Dnyr98FxtpI7JnS68zuMHvdgAREZED5YX81Pz4DBpao+Tf8L9EMvOJ+TLxxdq63Oa+rRP47bTlvTpO4ztPs/EnJ+3r9d8XG3+y/1FE7Uu/PWx9abDrzG7RFb+IiKSkvJCf+B3z+fzMEZQf+Zlu24atjwtXzOn1MfpT9HvafvPPPsFnh+7s1/6TQYVfRERSljGGH5w+mfriySz/8HWsnHN5t+2/svKYAUrWs6YPnnM7QqdU+EVEJKUVhNqfSsf8ISLB3G7bxmwG30pwGt9D9fUVv/CutX3azi0q/CIiktLGFmXz1tfbn6VHg7nUDZnUbfvKcIgf9qGn/9orS2hc8QLlN8319Nj+KvwiIpLyjh1byIPnz+LTRw1n/Yzze2xf3pLTp+Nsuf1jtG5aRqR2R6frd/z2kj7tN5Wo8IuIyKBw4dwxPHHxsXz/9MmEg3k9tr9n6xF9P1gs2uni2ld+3/d9pggVfhERGVS+fdpEtkye32O7JfWJT+RzqPXfKevztqlOhV9ERAaV/FCAxvzRCbWNOzB8ftOqF1l9eR4Vf/le/3eWAowdJJMK9MbcuXPt0qVL3Y4hIiJJ8tKG3Zz6y9fBWua8/KMu24UyYgDcf2TvBvcZCNMedK7+GmOWWWvn9txSV/wiIjIInTKxY5IdY1h1zMIu27XGfbTGfdy6cfIAJUt9KvwiIjKoNeeNYkvZGd22Wd3cc2fAdKHCLyIig9LrV5/IVz80DoyhatRxPbb/n12J9QvwupQv/MaYU40xrxhjfmWMOdXtPCIikhqOH1fEf3/uaF676sMJtf/HnmFJTjQ4JLXwG2MeMMZUGmM+OGT5GcaYNcaY9caY63vYjQUagRCwLVlZRURkcPrQ+CEJt/2JnvUn/Yr/D8BBD16MMT7gPuBMYBpwvjFmmjHmKGPMM4f8GQq8Yq09E7gOuDnJeUVEZBBads1JvHf8NTQUjOu23Ro968efzJ1ba182xow/ZPE8YL21thzAGPMIcI619jaguxEZaoDMrlYaYxYBiwDGjh3bj9QiIjLYzB5dSCQzn3Bmfo9tG6M+cv2xAUiVmtx4xj8K2HrA520dyzpljPmMMebXwB+Be7tqZ62931o711o7t7S01LGwIiLiLVesmUl5S7bbMVyT8p37rLVPWGsvt9aea6190e08IiKSmpZdcxJjJk5JoKXh4TTu4Z/UW/1d2A6MOeDz6I5lIiIifTZ7dCEv3X0dWVcMJyMeYdqyX3fZNm7NACZLLW5c8S8BJhljJhhjgsB5wNMu5BAREY8xxpBbMpzWrGJiGYEu25W3ZrO0vnAAk6WOZL/O9zDwOjDFGLPNGLPQWhsFvgY8C6wCHrPWrnDoeAuMMffX1dU5sTsRERmEnrv8eKwvwPKTup5UJ2YzuHvrxAFMlTo0SY+IiHjOz1/awLVPr8QfbmTm63d22a44EOa/Jr8/gMn20yQ9IiIiDvnmKe1X89FgLu9+6Nou2+2JBAcqUspQ4RcREU+LBnO7XR/33o3vbqnwi4hIWvvRxqluRxhQnir86twnIiKHGlUQ4u2TbqCm5MhO15e35HDhijkDnMo9nir81trF1tpFBQUFbkcREZEUseF7H8Fm+Cmf/gW3o6QETxV+ERGRvSaV5ACQ6ffx/OXHu5wmdajwi4iIJ73/7VNouu1MAD42uX0Ol1WzL+uy/aUrZw1ILrep8IuIiCdl+n1kBw8emb45byQbpn2+0/Zh6xuIWK5T4RcRkbTwlwv3duDrepz+OzaXDUwYF3mq8KtXv4iIdOVzM0dyysRiWnKGdtnmvcYCfrCh897/XuGpwq9e/SIi0p3JpTm0ZRezcs5Xu2yzuTV7ABMNPE8VfhERke4cO6Z9Rr6W3GEuJ3GPCr+IiKSNS48by9rrT+P4cUVuR3GNCr+IiKQNYwyTSnN5/eoTqTjjRrfjuEKFX0RE0lL5zfPdjuAKFX4REUlLAV96lkBP/a31Op+IiDjhTztHux0haTxV+PU6n4iI9Mbaoy5gx7hTDlv+z+qu3/Uf7DxV+EVERHqj/vffYOe4U6gacfC0vNalPANBhV9ERNLaESU51JQcPFqf7WZY38HO33MTkfS2detWmpub3Y4hHpOdnc2YMWPcjiHAmutO4403Qnzj639yO8qA0BW/iIikNb8vgxnTpx+2/Kmq4bxd770+Y7riF+mBrspEvK+wsJBwZj7Btvp9y/5f5SgAHpq+zK1YSaErfhEREWBL2ZluRxgQnir8eo9fRET66gvzP+52hAHhqcKv9/hFRKSvrJff4TuAnvGLiAywcMV64q2NbseQQ5TWrWFsqOs3eHaHA5QEI44dr3Xzcsf21Rsq/CIiAyy8ezOxht1ux5BDLBhaT22wtcv1ARN39HhtO9c4ur9EebLwN0diLN+u5/zinNxMP2UlOW7HEJEkmjYsv9Plj+wcSRzD/JJdA5woOTxZ+Otbozy7utLtGOIhI/JDKvzimGDJOOK5xW7HkE7sCoe6XNZc6ePSUVscO1bmiCmO7as3PFn4M30ZjB+S7XYM8YhNq95l99Zmnq1e6XYU8ZCioiLmzZvndgw5xJbWrmvHltZsRwt/aNwsx/bVG54s/EXZAc49ZpTbMcQjntn+DlVVDezY0eB2FPGIvOoVBLMzqItvdDuKHOJb86fz8suvDMix6t54dECOcyhPFv6ttS18/cn33Y4hHlGzsZF8fCyYPsztKOIR2a0ryKGZ8O5NbkeRQxw9tpiVgbYBOZZb59+ThT8as1Q1DsyJE+8LB/LJDRnGjx/vdhTxCJv/cUK+OHnDhrsdRQ7RUlnJBw/8o5sWOxw7Vt5Rn3BsX3B9wi09VfiNMQuABWMmlPGdj0xyO454RMXObHyxsNsxxENs3XbisWba4nr7KNUUADkzT4Gda8ir3ZTUY+l1PgdYaxcDi6fNPOYyt7OId7z9+qs01ewmP+Spbxdx0YTKFxkSjJA1bZrbUaQTZ5XWs6Gqkuys5A6y1Lz+taTuvyue/Emm1/nESW+8/S6xxhpKczPdjiIeURLfQigzRrSh1O0o0ok5RTF2ZvnJa40S82fiiybn0XG0oSop++2JJwu/iJNyikohGGD4kCy3o4hH+OtL8YVi+PNU+FNRHnD+6cO566UNYAxDt7+5b92KhlzGhprJC/R/FD+3zr8nC38kFmdHfdfDLor0Ru70EynJhIuOHeN2FPEIW30aIV+cInXuS2lPLXkJgCPrNhy27scTV/V7/0UnfaXf+9jv7oRberLwi4ikMjNkHL7sbEJj9MtkKlvwsSA/e3EDOd0M6tMfGsDHQcU5Qb4yb6zbMcQjVr67m2hzA5s2bXI7iniEv2EHBVkBSuMz3I4i3Tg+UMHUaHmnM/bFLPhM//av2fkclB3wMWtUgdsxxCOK42U0N3c9VadIb8U37iAzVufa61ySmI8UWBZOaKX8/cMfHVsH9q/X+URS1M6dO6mpqXE7hniIvxEKsgoY6dIkLZK4GUeV8Nrjh0/cs9fGliwmZLX0ad9uTdKT4cpRRUREBoGuruw3tbQ/91/R2PlUvqlMV/wiPRgxYgQFBXp0JM7Rrf7BI1Bdw/Dg4bf6X6gu4a26GCMz+3a1D7rV76hdDW389IV1bscQj2jcvYvCgGW+JukRh5iCUWT44mTqdb6UF2naza5w17f6c/3RPu/brVv9niz8Ik7KLRlOaX6IKVP06pU4Rc/2B4u2+p0sKT6NYQcM4nOgqDV8srSiT/vW63wO2DtJT1lZGdd9VJP0iEhqalzxArGG3W7HkATkbq3h5MA6CvKru2xz79YJfG3Mxl7vu+6NR/sTrc88Vfj3TtIzcdrMyx59Z7vbESdTLgAAGxJJREFUccRDSnKCfHSyhlcVZ7RseptozTa3Y0gCQrubONK3KykT9miSHhGRNJE1fjaxYg0yNhhkFzWx6tUtFLT0VPh7P+FOdtkJfQvVqTQfsrcoO8C5x4xyO4aISKdyp3/U7QiSoNnA6CU7eeGZrm/1t+v9rf6C48/tU6bOnZdwS08WfhEnbd26VSP3ieOys7MZo7H6B4W9U3KHg3kEww0up+k/DeAjIiLSjZkzZwKwZ9hMl5M4Q1f8Ij3QVZlIejv99NM55phjKL3132BgxJZX3Y7UL7riFxER6UFJSQl1t57FjgmDv3+GCr+IiEgC8kMBALaUnelykv7RrX4RkQEWrlhPvNX598Il+aZGy2FYCWO39b/Db+vm5Q4k6j0VfpEeqFe/OM1WbybkizNcY/WLC3SrX0REJI3oil+kB+rVL04LV/h0q3+QuvWsqXzv76vdjtEvuuIXERFJ0JlTB//03Cr8IiIiacSTt/qbIzGWb69zO4Z4RMXO7fhiYcYUZrkdRTxCnfvETbriFxER6aV10893O0KfefKKPzvgY9aoArdjiFfo/5I4rHHFNmINNbTt1J3JwWh8dBv4YXiwdd+y5/eUMDbUzJScxF/9bdu5JhnxeuSpwm+MWQAsKCsrczuKeIje4xenxVYuIdhWR0FhodtRpA9GxioBKAm07Vu2JxJkTyTYq8If3r3J6WgJ8VTht9YuBhbPnTv3MreziIh0JWP0bPy+OHl6xj8ovfZ8CIAjm5Z0snZHwvvJO+oTDiUCuD7hlp4q/CLJoPf4xXlT3A4g/bDavxWAnNbsfu0nNG6WE3F6TYVfRGSANa54gVjDbrdjSB+d0boUgAn51f3aT90bjzoRp9c8Wfj1Op84aeW7bxNtbmBYXqbbUcQjssufJ4dW3U0apDp7xt8XesbvoLZInDWVGg5TnLG1toVAOKzCL46JZZdgfTGCJePdjiJ9sMPXfrcmM9K/nwlunX9PFv7MQAZThua6HUM8YvcHTbQ21bBjR4PbUcQjhm1eRtDXSjOVbkeRPpgZKQcgY0IZ2Q07yW6q6NN+mte/5mSshHmy8Is4qb62hqaa3cSa9O0izshtbCIzGHE7hvTTx6aN5p8rOajwb23JZExW/x4BJJsnf5LpVr84qWjybEabKPPGFbkdRTwi9l6cYFsdfr3HPyhVZdQCcNSUKbT683hj24p96xpjfiCxwu8vHJmMeD0f15WjJlljOMprm/rX21Jkr5bqWkoyUeEX52xdRrSpiqYszf8wGM2NVBDyZ9C0upL83fVMzd7/GHCIP/Gr/abVLyYhXc88WfhFRESS5fRJJZiOmW58GcbdMH3gycIf8GUwMj/kdgzxiHdfeZXGplr+sEJXZ+KMWbVbKA5GmDLVnVu90j++A74uyMilKbZy3+eITXzuO1+OO3cRPVn4I7E4O+pbe24okoCK7ZuxTXXE64JuRxGPOIFyQplhWjY2uR1F+isaY0Jo/3nM8yfeabNl47JkJOqRJwu/rvjFSfWTZ5DRUs+UoTluRxGPKFnzT0I2Sqxxj9tRpJ+MhcwsP/5w+8WmH5vwtm6df08Wfl3xi5O2b1oPTXXE92gAH3HGh+ItxIlgYwG3o4gDfFgyTHvB780TfxuLJidQDzxZ+JvCMZZuqXU7hnhE2+49+MIN5FjdRRKHZEYhMw7xmNtJxAE+4mTQ+8Lv1vn3ZOGPxuLsbkrtARRk8AgUjCY72sDUSSVuRxGPaKtcStC04i8a4nYUcUBVRT2haPsr5Dn+KJnEE9rOX+Rk5861Cbf0ZOHP9Ps4oljPY8UZtdF6/G317Nrlzm058Z5wvI1ohv4/iTs8WfhHFYa4bf40t2OIR/zXsiANNdDS0uJ2FPGIUvYQzIwTrdD8D14wmjjGn3invr2iFRuSkKZnniz89a1Rnl2tyS/EGWtXrSQWbmUQjtMhKSo63uKPAJmJv/MtqSseB0P78/pendEMJ89/4v0FPFn4q5vCPPz2NrdjiFe0NOMnnuBTO5GeZWDbO4HZ3l8lSuoxe89nb7l0/j1Z+MOxOJtqdFtWnDEGsIBevBKnBOno/W31nN8L+nzd7tL592Tht0Asrt+kxSkZQBxjdK9fRAY/TxZ+EUf5/BCL4/PpeayIDH6eLPxxC80RPZEVh8TCQJxoVP+nRORwASDxEfrd58nCbwBdnIlzVPDFWWEgE8B48kdwGrLYjl71vXog6Oj5T7y/gCf/12UYyM/05F9NRDygLhYgz8TIzC5wO4o4wMbjtLTVARA08YQ7AhtHz3/iE/6kfHU0xmQAtwD5wFJr7YM9bROzUNeq3rLijBFuBxDPiVpDzBqMT++KeIE1MWK2/Vo/3otLfrfOf1ILvzHmAWA+UGmtnXHA8jOAuwAf8Ftr7X92s5tzgNG0/zqT8Mv5UXXqF5EUtaE1l+JghNHjZrkdRRwQCYfZsmMJAMOCreQEEnviH3L0/P8j4ZbJvuL/A3Av8NDeBcYYH3AfcDrthXyJMeZp2n8JuO2Q7S8BpgCvWWt/bYx5HHghyZlFDpYRhHhEvfrFMS3xIG348OUUuR1FHGCDEZpi7eU0YhP/OeHW+U9q4bfWvmyMGX/I4nnAemttOYAx5hHgHGvtbbTfHTiIMWYb7X1hoJteVsaYRcAiAIZN7G90kX3C2cUEwo2UFGW7HUU8YkjWDoqzINZU43YUcYABcnztj5cDJvHOwG6dfzee8Y8Cth7weRtwXDftnwDuMcacBLzUVSNr7f3A/QBm+CTd6BfHNBaOJ9BWz9FHOzmFpqSzvNrdZAcjuuL3EF3xO8ha2wwsdDuHpC9ftIVAtIWGBs2kJs4Itu3BHw0T2b3Z7SjikJJAGwBZGYlPluPW+Xej8G+nffjzvUZ3LHOMP8NQkqvesuKMLL+PkE3535FlEMmKN5NJG9GG3W5HEYccMaKI7Q1tBCOtCW/j1vl346fZEmCSMWYC7QX/POCLTh4gP9PPaWUlTu5S0tjOjSH8reGeG4okqNI/DBuMMXrCjJ4by6CQBfxj2TZmtrxFQcfVf4/bTJjjYIK1CbdM9ut8DwOnAiUdnfRustb+zhjzNeBZ2nvyP2CtXeHQ8RYACwpGl1Gam+nELkVoyM/HH7Dk5eW4HUU8oiYwlUB+kPxjznI7ijjoxQ+Wktm8irFZbexqDTI81P0FQ/4x5zh49IcTbpnsXv3nd7H878Dfk3C8xcDi4gnTLqtqTOw3LpGetISjBDROvzioyV9AcyiPzBFT3I4iDtrk38XLtSWUt+SyrS3Ejyeu7ra9W+ffkw8u22JxNle3uB1DPCJcV09mpB717ROnmHgjuZkaF8KbDNvastwO0S1PFn4RJ2XklRKIZzF27FC3o4hHZLdlM7Qo1+0YkqY8Wfhzgj7mji10O4Z4xNbtWQTbNPeDOCcrXEOouYG2nWvcjiIOGh/dxvBg4r363Tr/nir8ezv3DRtXxgnjh7gdRzzigx0jobme0hKN3CfOyPNVkpetW/3iDk8V/r2d+wrGHXnZH97a4nYc8YhYUwFDQ4VccMJkt6OIR9jqUYR8cTKHDXc7ijhok38XWeFQwu3Vuc9BkVicnfWJ324R6VZjGzbai7k2RURSmCcL/5CcIOfPHu12DPGIinU1ZNs2Nm3a5HYU8Qh/ww4KsgIM1xW/uMCThV/ESTlFpRQGLOPHD3M7iniErTaEfBobQtzhycIficXZoVv94pCWxjZiGghSRDzCU4V/b6/+MRPK+Mq8sW7HEY+o2OnDF9NY/SLiDZ4q/Ht79c+dO/eyWaMK3I4jXqH/S+I4DdXrRav9W8lpTfy139C4WUlM0zVPFf69miMxlm+vczuGeEhupp+yEk3SIyKDnycLf1skzprKRrdjiIeU5ARV+EWkRxunnMOENU+5HaNbniz8mYEMpgzVONjinNxMT36riIjDqofPUuF3Q3bAh57xi4iIm5bVFzAnP/UeO3uy8Nc0R3j0ne1uxxAPKckJ8tHJpW7HEJFBpCEWcDtCpzw1S4QxZoEx5v6m5ma3o4iISNqzbgfolKeu+A98ne/cY0a5HUdERCTleOqKX0REJFWk6tReKvwiIiJpxFO3+kWSYevWrTSr34g4LDs7mzFjxrgdQ9KQCr+IyACz1ZuJ1cVpje9xO4o4aGq0HICxofYLhVxftNv2rZuXJz1TZ1T4RXqgqzJxWrjCR7xVo4uKO1T4RUQGWHBYmdsRJAlW+7cC7JuopynWfYnVJD0OOHBaXk3SI07SJD0i0lsrmvI5uSj1Hud4qle/tXaxtXZRbp7G6RcREXetbU7NWuSpwi8iIiLdU+EXERFJCg3ZO2A0O5+IiEjnPFn4RZykAXwkGTSAj7hFt/pFRESSoLGH1/nckpqp+qk5EtPrfOKcjHxyi4fodT5xTLhiPfHWPbRuTr1XvaTvDh25rycauc9BbZE4ayo1KpY4Y9OqdwmEm5k+Is/tKOIR+e/8lqxIHcUlJW5HEQdd3LQNgCFDdiXUvurvtyczTpe8WfhjcTZV65msOKOivpVsG3E7hnhIPDMf68/An1fqdhRxUE1GAwCBBG/xu3X+PVX4947cN2xcGeOHZLsdRzxi/HHHUZIT5KOT9UNanBGeNVFj9XvQU0teAuDIug37lp1RUtVl+6KTvuLg0e9OuKWnCr+1djGweO7cuZede8wot+OIR7T36q9mzZpqt6OIR8Q3vkZmrJniYt3q95Jrp8V5f1c9kWBrQu3bdq5JcqLOearwiyTDzp07qampcTuGeIi/EQqyChg5YorbUcRBF4yYwp/f3sbfnn80ofaZLp1/Txb+muYIj76z3e0Y4hEfLFkDzfUcUaLHR+Kc0lAhc12anU2Sp3FzPlta9/+suHvLEVw9trzTtpqdz0Hq3CdO2t3YRiASdTuGiAwCxhz8eWlDkTtBuuHJwi/ipJLxUygMWE6YPsztKOIh2dm6g+RFNjWH5z+IJwt/pi9DvfrFOUOOoCQnyBT16heHNK54gdj23dTpiaTnDN1Qwbz8xDoC172RWF8Ap3my8BdlB1CvfnHKW2+9Rc3GGp7d6HYS8Yrs8ufJoVVj9XtQqGEPJYG2hNqGd29KbpgueLLwq3OfOEmd+8RpQ2pqIKh+I+IOTxZ+de4TJ6lznzgtEK7HF2slWrvD7SjisMzmWgr9iY306db592Th1zN+cVJjbiY0J3brTkQk1Xmy8IuIpLLmvHFkZWeQXTbH7SjisLpIBetbXjlkaefD9maXneDgkdN0yN4Dx+oXcUpuUTGBnCxGanY+cYg/D/KyAq6N3CbJ07wzh13hUEJtNXKfA/aO1T9x2szL3M4i3jH+yJmapEccFa5Yr0l6xDWeKvx7ZQYymDI01+0Y4hEVO7dja8KapEccY6s3E/LFGT5suNtRJA1luB1AREREBo4nr/hFRFKZGTIOX3Y2IQ3g4zkNh0zS0x1N0uOg7ICPWaMK3I4hHvHW9jXU1NawqdbtJOIlRUVFGrnPg+wgGKzfk4VfxEkjRoygoEC/SIqzNEmPuMWThb85EmP59jq3Y4hXZOSTWzyEspIct5OIR7T36t9D6+Y9bkcRh+Xt2crYUGIjx7ZuXp7kNJ1T5z4REZE04skrfj3jF5FUFhymQca8qmFTXsp37tMVv4iISBpR4RcREXFI6vfpV+EXERFJKyr8IiIiaUSFX0REJI2o8IuIiKQRT77OpwF8xGm5mX4N4CMinqArfhEREYcMgqH6vXXFb4xZACwoKyvTAD4iIiKd8NQVv7V2sbV2kSZUERER6ZynCr+IiIh0T4VfREQkjajwi4iIpBEVfhEREYfkh1K/z7wKv4iIiEMmFqf+eB8q/CIiImlEhV9ERCSNqPCLiIikERV+ERGRNKLCLyIikkZU+EVERBy0avZlbkfolgq/iIiIg5rzRrodoVsq/CIiImlEhV9ERCSNqPCLiIikERV+ERGRNKLCLyIikkZU+EVERNKICr+IiEgaUeEXERFJIyr8IiIiaUSFX0REJI2o8IuIiKQRFX4REZE0osIvIiKSRvxuB+iJMeYk4Eu0Z51mrT3B5UgiIiKDVlKv+I0xDxhjKo0xHxyy/AxjzBpjzHpjzPXd7cNa+4q19qvAM8CDycwrIiLidcm+4v8DcC/w0N4FxhgfcB9wOrANWGKMeRrwAbcdsv0l1trKjq+/CCxMcl4RERFPS2rht9a+bIwZf8jiecB6a205gDHmEeAca+1twPzO9mOMGQvUWWsbkhhXRETE89zo3DcK2HrA520dy7qzEPh9dw2MMYuMMUuNMUurqqr6GVFERMSbBkWvfmvtTdba13poc7+1dq61dm5paelARRMRERlU3Cj824ExB3we3bFMREREksyNwr8EmGSMmWCMCQLnAU+7kENERCTtJPt1voeB14EpxphtxpiF1too8DXgWWAV8Ji1dkUyc4iIiEi7ZPfqP7+L5X8H/u708YwxC4AFZWVlTu9aRETEEwZF575EWWsXW2sXFRQUuB1FREQkJXmq8IuIiEj3VPhFRETSiAq/iIhIGvFU4TfGLDDG3F9XV+d2FBERkZTkqcKvzn0iIiLd81ThFxERke6p8IuIiKQRFX4REZE0osIvIiKSRjxV+NWrX0REpHueKvzq1S8iItI9TxV+ERER6Z4Kv4iISBpR4RcREUkjKvwiIiJpRIVfREQkjXiq8Ot1PhERke55qvDrdT4REZHuearwi4iISPdU+EVERNKICr+IiEgaUeEXERFJIyr8IiIiaUSFX0REJI14qvDrPX4REZHuearw6z1+ERGR7nmq8IuIiEj3VPhFRETSiAq/iIhIGlHhFxERSSMq/CIiImlEhV9ERCSNqPCLiIikERV+ERGRNOKpwq+R+0RERLrnqcKvkftERES656nCLyIiIt1T4RcREUkjKvwiIiJpRIVfREQkjajwi4iIpBFjrXU7g+OMMVXAZgd3WQA48Y5gX/aT6DY9tetufWfrump/6PISYHcC+ZIh1c9Lf9r0Znlny9LxvLj1vdLVci+el8H0M6yr5al0Xpz8Xim01pYm1Npaqz89/AHud2s/iW7TU7vu1ne2rqv2hy4Hluq8ON+mN8u7WJZ258Wt75V0Oi+D6WfYYDgvbn2v6FZ/Yha7uJ9Et+mpXXfrO1vXVXun/i2ckOrnpT9terM8lc4JuHde3Ppe6Wq5F8/LYPoZ1tXyVDovrnyvePJWvwwcY8xSa+1ct3PIwXReUpPOS2pKt/OiK37pr/vdDiCd0nlJTTovqSmtzouu+EVERNKIrvhFRETSiAq/iIhIGlHhFxERSSN+twOItxhjTgVuAVYAj1hrX3Q1kABgjMmg/bzk0/7O8oMuRxLAGHMS8CXafxZPs9ae4HKktGeMGQvcDVQDa621/+lyJMfpil96ZIx5wBhTaYz54JDlZxhj1hhj1htjru9YbIFGIARsG+is6aSX5+UcYDQQQeclqXpzXqy1r1hrvwo8A+iXsSTp5ffKUcDj1tpLgGMGPOwAUK9+6ZEx5mTai/lD1toZHct8wFrgdNoLyRLgfGC1tTZujBkG/Nxa+yWXYnteL8/LJ4Eaa+2vjTGPW2s/51Jsz+vNebHWruxY/xiw0Frb4E5qb+vl90oF8DjtFzF/tNb+3pXQSaQrfumRtfZl2m97HWgesN5aW26tDQOPAOdYa+Md62uAzAGMmXZ6c15o/8FW09EmjiRNL8/L3lvLdSr6ydPLc/IV4CZr7UeAswc26cBQ4Ze+GgVsPeDzNmCUMeYzxphfA38E7nUlWXrr9LwATwCfMMbcA7zkRrA019V5AVgIeO6qchDo6pz8A7jaGPMrYJMLuZJOnfvEUdbaJ2gvMpJCrLXNtBcYSTHW2pvcziD7WWs/ADz9KExX/NJX24ExB3we3bFM3KXzkpp0XlJP2p4TFX7pqyXAJGPMBGNMEDgPeNrlTKLzkqp0XlJP2p4TFX7pkTHmYeB1YIoxZpsxZqG1Ngp8DXgWWAU8Zq1d4WbOdKPzkpp0XlKPzsnB9DqfiIhIGtEVv4iISBpR4RcREUkjKvwiIiJpRIVfREQkjajwi4iIpBEVfhERkTSiwi8i/WaMKTTGXHHA55HGmMeTcJwXjTFznd6vSDpR4ReRhBhjupvboxDYV/ittTs09a9IalLhFxnEjDE5xpi/GWPeNcZ8YIw5t2P5HGPMS8aYZcaYZ40xIzqWv2iMucsYs7yj/byO5fOMMa8bY94xxrxmjJnSsfxiY8zTxpj/A14wxuQaY14wxrxtjHnfGHNOR5T/BCZ27PdnxpjxxpgPOvYRMsb8vqP9O8aY0w7Y9xPGmH8YY9YZY25P8K/95U7y/9AY88eOv8M6Y8xlTv0bi3iNZucTGdzOAHZYa88GMMYUGGMCwD3AOdbaqo5fBn4CXNKxTba1dpYx5mTgAWAGsBo4yVobNcZ8DLgV+GxH+9nA0dba6o6r/k9ba+uNMSXAG8aYp4HrgRnW2lkdOcYfkPFKwFprjzLGTAWeM8ZM7lg3CzgGaAPWGGPusdYeOFVqZzrLD3A0cDyQA7xjjPmbtXZHov+QIulChV9kcHsfuNMY81PgGWvtK8aYGbQXw+eNMQA+YOcB2zwMYK192RiTb4wpBPKAB40xkwALBA5o/7y1trrjawPc2lF047TPXz6sh4wn0v6LCNba1caYzcDewv+CtbYOwBizEhjHwXOkd6az/ABPWWtbgBZjzL+AecBfe9iXSNpR4RcZxKy1a40xs4GzgB8bY14AngRWWGs/1NVmnXy+BfiXtfbTHVfrLx6wvumAr78ElAJzrLURY8wmINSPv0LbAV/HSOxnUmf5u1suIgfQM36RQcwYMxJottb+CfgZ7bfl1wClxpgPdbQJGGOmH7DZ3n4AJwJ1HVfcBeyfi/zibg5ZAFR2FP3TaL9CB2ig/a5BZ16h/RcGOm7xj+3I2N3f66G9z+870Vl+gHM6+hMUA6fSPu2qiBxCV/wig9tRwM+MMXEgAvyHtTZsjPkccLcxpoD27/NfAHunHG01xrxD++38vc/9b6f9Vv/3gb91c7w/A4uNMe8DS2nvG4C1do8x5t8dHfr+F7jvgG1+Cfx3xzZR4GJrbVvHY4iuHA109Xy+s/wA7wH/AkqAW/R8X6RzmpZXJI0YY14EvmWtXep2lq4YY/KB31lrP9+LbX4INFpr70haMBGP0BW/iKQUa209kHDRF5He0RW/iIhIGlHnPhERkTSiwi8iIpJGVPhFRETSiAq/iIhIGlHhFxERSSMq/CIiImnk/wNTyMkHDlBp/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f3f0cb67da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gs = plt.GridSpec(nrows=2, ncols=1, height_ratios=[10, 2])\n",
    "plt.figure(figsize=(8, 15))\n",
    "\n",
    "SMOOTH = False\n",
    "\n",
    "ax1 = plt.subplot(gs[0])\n",
    "for i, cond in enumerate(conditions):\n",
    "    \n",
    "    ref_point = 200000 // binsize\n",
    "    norm = 1 #scalings[cond][ref_point]\n",
    "    \n",
    "    # cis P(s)\n",
    "    x = np.arange(0, len(scalings[cond]) * binsize, binsize)\n",
    "    y = scalings[cond] / norm\n",
    "    if SMOOTH:\n",
    "        x, y = coarsen_geometric(sums[cond], n_valid[cond], 100)\n",
    "        x *= binsize\n",
    "        x = x[:-1]\n",
    "    \n",
    "    plt.plot(x, y, \n",
    "             color=colors[cond],\n",
    "             label=long_names[cond])\n",
    "    \n",
    "    # average trans levels\n",
    "    for _, row in trs_exp[cond].iterrows():\n",
    "        plt.axhline(\n",
    "            (row['balanced.sum']/row['n_valid']) / norm, \n",
    "            xmin=i/len(conditions),\n",
    "            xmax=(i+1)/len(conditions),\n",
    "            c=colors[cond], \n",
    "            alpha=0.25)\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.ylabel('P(s)')\n",
    "plt.xlabel('separation, bp')\n",
    "plt.legend()\n",
    "plt.gca().set_aspect(1)\n",
    "xlim = plt.xlim()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
